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INTRODUCTION 

Tensile testing of nominally prismatic round and flat bars 
has long provided a basis for evaluation of the properties of 
metals. The primary data from such tests, load, extension and 
minimum deformed cross-sectional area may be used to infer stress- 
strain relations taken to be characteristic of the material tested. 
Inference of such material property relations is appropriate so long 
as the deformation and associated stress field may be assumed 
uniform over the cross-section in some region or gage length of 
the tensile bar. 

For ductile metals the assumption of uniformity of the 
deformation and stress field clearly breaks down for extensions in 
excess of some material dependent critical value. Once the 
extension exceeds this value the axial load decreases and subsequent 
prescribed overall specimen extension is observed to concentrate 
in a highly strained local neck. Within this necked region 
inhomogeneous deformation and a complex multi -axial stress field 
must. exist. The phenomenon is referred to as tensile instability. 

Evaluation of the validity of stress-strain relations inferred 
from post-instability tensile test data as well as study of tensile 
failure by fracture and/or rupture requires detailed knowledge of 
the stress and strain fields existing in a necking tensile bar. 

This need has motivated extensive experimental investigations, most 



notably by Bridgeman [1 ] and Nadai [ 2], as well as a variety of 
approximate analyses incorporating idealized models of material 
behavior. Typically these analyses neglect either elastic deforma- 
tion or work hardening, or both. 

The absence of a complete theory suitable for analysis of 
problems of general finite deformation of elasto-plastic conti nua, 
such as necking in metal tensile bars, has provided the motivation 
for development of such a formulation. 

Analysis of finite deformation of elasto-plastic materials 
requires explicit consideration of nonlinear effects arising from 
both inelastic material behavior and deformation magnitude. While 
the recent literature contains many examples of analyses incorporating 
either material or geometric nonlinearity, few attempts have been 
made at solution of the combined problem.* Excellent surveys of 
these efforts are provided by Marcal [4] and Stricklin et. al . [5]. 
These previous developments typically employ finite element 
techniques based upon energy principles and are primarily intended 
for analysis of problems involving large deflections of plates 
and shells in which local strains are small. The application of - 
these analyses to the tensile instability problem is of questionable 

*0den L3J provides extensive development of theory and solution 
techniques for finite deformation of hyperelastic materials. The 
approaches employed, while incorporating both material and geometric 
nonlinearity, do not admit application to the elasto-plastic case. 
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value since the necking process involves large displacements, strains 
and rotations distributed over a continuum of arbitrary shape. •The 
kinematic assumptions underlying analyses of plates and shells are 
of limited validity under these conditions. Solution of a general 
elasto-plastic continuum problem requires a formulation appropriate 
for analysis of deformations of any magnitude irrespective of the 
configuration of the deforming solid. 

The adoption of what is herein termed the. note vtewpoint* 
toward the mechanics of finite deformation of an elasto-plastic 
solid has led to the development of a complete theoretical 
formulation of the problem. Rather than seeking equations 
governing the total deformation attention is restricted to the 
time rates of the independent variables, stress rate and velocity. 
Equations are derived governing the time dependent velocity field 
in a deforming elasto-plastic solid. 

The formulation differs from previous developments in two 
fundamental respects. 

*A similar viewpoint is taken by Cowper and Onat [6] in establishing 
admissible solutions for tensile necking in plane strain but they 
do not attempt full solution for deformation and stress histories. 
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1. The entire development proceeds in an Eulerian or spatial 
reference frame rather than the Lagrangian, or material, 
frame usually employed in nonlinear analyses of plates 
and shells. 

2. Constitutive equations for finite elasto-plastic deformation 
are obtained by generalizing those of the infinitesimal 
theory in the spatial rather than material frame. The 
fundamental features of the infinitesimal theory are 
preserved without introducing problem dependent deviations 
associated with the use of material frame stress tensors. 

A complete initial- and boundary-value problem is posed in which 
finite elasto-plastic deformation is viewed as a time dependent 
process. The formulation reduces to well established results in 
the limit of infinitesimal deformation. 

The governing equations of the finite problem are distinguished 
by their quasi -linear nature. This feature, which follows directly 
from adoption of the rate viewpoint in an Eulerian frame, enables 
the use of an incremental technique for accurate and efficient 
numerical solution of finite deformation problems. Finite element 
solution capability may be developed directly from the governing 
differential equations. 

Numerical procedures have been developed for analysis of finite 
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deformation under conditions of plane stress or plane strain. The 
capabilities of these procedures have been investigated by considering 
a number of problems of homogeneous finite deformation for which 
analytic solutions are available. Comparison of numerical and 
analytic results for these problems indicates that accurate numerical 
solutions can be obtained for problems involving dimensional changes 
of an order of magnitude and rotations of forty-five degrees. 

The numerical analysis has been employed in an investigation of 
symmetric necking in flat tensile bars of elasto-plastic material. 
Solutions are obtained for the limiting cases of plane stress and 
plane strain extension of bars containing a small initial geometric 
imperfection. Full histories of neck geometry and internal stress 
and deformation fields are obtained. The development of inhomogeneous 
internal fields as the necking process proceeds is clearly demonstrated, 
as is elastic unloading of previously yielded material in regions 
outside the neck. The solutions also provide a vehicle for assessment 
of the validity of stress-strain relations inferred from tensile 
data over the full range of a test from initial yield through the 
development of a significant neck. 
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I. THEORETICAL DEVELOPMENT 


The Rate. I /tewpotnt 

The equations describing finite deformation of elasto-pl asti c 
solids may be derived in what is termed a rate form. That is, atten- 
tion is focused not upon field quantities such as stress and strain 
but rather upon their rates of change with respect to time. The 
approach is conceptually analogous to that employed by Swedlow [7 ] 
for infinitesimal deformation. Even with this analogy, however, 
analyses of infinitesimal and finite deformation are operationally 
distinct. 

In analysis of infinitesimal deformation stress and strain 
tensors as well as all governing equations are referred to a single 
configuration of the body. Either deformed or undeformed states 
may be employed as they are by assumption indistinguishable from 
one another. Thus time derivatives of field quantities reflect only 
changes in component magnitudes with respect to an invariant frame 
of reference. Should the deformation be regarded as finite, however, 
deformed and undeformed configurations must be distinguished. Time 
derivatives of field quantities such as stress and strain must 
reflect changes in the fundamental reference frame provided by the 
deforming configuration of the body. 

The present theory incorporates equations of elasto-plastic 
behavior and stress equilibrium which reflect the foregoing implica- 
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tions of finite deformation. The equations are assembled to define 
a complete initial -and boundary- value problem for the time dependent 
velocity field in a deforming body. 

The development of the constitutive and equilibrium equations is 
predicated upon the character of certain tensorial measures of stress 
and strain and their time rates of change. In the next few sections 
these quantities are defined and discussed. Subsequently the field 
equations are derived, boundary and initial conditions developed 
and the full velocity problem is assembled. 
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1 . 1 Fundamental Concepts 


General concepts of nonlinear continuum theory are developed 
below to the extent necessary to support the ensuing analysis. More 
extensive discussions* may be found in Eringen [8 ] and Truesdell 
and Toupin [ 9 ]. 

Descriptions of continuum deformation and loading are developed 
in a fixed or laboratory reference frame. Deformation, deformation 
rate, stress and stress rate tensors are defined. All definitions 
are subject to the constraint of material objectivity, or spatial 
invariance, which requires that the analysis be independent of 
rigid motion of a deforming continuum. 

Continuum Motion : Consider a three dimensional body whose 

undeformed reference configuration is B 0 with boundary 3B 0 . The 
deformed configuration of the body is B with boundary aB. Material 
points in the reference state are located by coordinates X* in an 
orthogonal curvilinear coordinate system. In the deformed state 
these points are located by coordinates x 1 in the same system. 

Thus the time dependent deformation may be given as the mapping** 

*Much of the discussion in this section has been abstracted from 
the texts cited. 

**General tensor notation is employed. Repeated indices in subscript- 
'Superscript pairs imply summation . over 1 ,2,3. A comma denotes 
partial and a semi -col on covariant differentiation. 
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from B 0 to B in Figure (1). 


x 1 ' = x i (X K , t) (1) 

The nature of the motion (1) is limited only by the constraint of 
material continuity: 

, 0 = |x; K 1 > 0 (2) 

In the special case where the mapping of B 0 to B given by (1) 
is such that angle and distance are everywhere preserved the motion 
is termed rigid. In this case (1) admits the representation 

x 1 = Q 1 |< X K + b 1 (3) 

In (3) the tensor Q 1 ^ is orthogonal, (4)*, and both and b 1 are 
functions of time.** 

o^W'j (4 > 

Reference Frames and Time Derivati ves : We may def i ne a 

i K 

velocity field v for the motion (1) by noting that the X are 

constant and differentiating (1) with respect to time. 

*In (4) overscript T denotes transpose. Thq mixed metric tensor 

g 1 . is equivalent to the Kronecker delta 6 1 .. 

J J 

K l 

**In cartesian coordinates Q . is a rotation and b a translation. 

In general coordinates no simple physical interpretation is possible. 
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(5) 


v 1 = axVat = v i [x i (X K ,t),t] = v n (X K ,t) 

Note that the coordinate dependence of v 1 may be given in terms 

"i K 
of either x or X . 

The choice of coordinate dependence in (5) is characteristic 
of all tensor field quantities employed in finite deformation 

V 

analysis. Analysis employing the X is termed material or Lagrangian 
while that employing the x 1 is termed spatial or Eulerian. The 
present theory is developed in Eulerian form. 

In differentiating Eulerian field variables with respect to 
time, the time dependence of the x 1 must be fully taken into account. 
Thus if <j> ( x 1 , t ) is a tensor field of any order its time derivative 
<(> is 

^(x 1 ,t) = 3<t>/3t + 4> . . v 1 (6) 

» I 

The derivative (6) is the total of material derivative of <j>. 
Unless specific exception is noted it is the only type of time 
derivative employed in the present theoretical development. 

Material Objectivity : Any theory attempting to describe physical 

phenomena must be independent of the observer of an event. This 
constraint is known as objectivity. A familiar example of its 
significance is the restriction of Newtonian mechanics to inertial 
reference frames. 
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In continuum analysis objectivity is most conveniently expressed 
as the requirement that any theory be independent of rigid motion. 

The constraint may be expressed operationally by considering tensor 
fields referred to deforming bodies whose motions differ by a rigid 
component. Objectivity requires that components of such tensors 
be related by the transformation relating the two motions. 

Consider, for example, the velocity fields corresponding to two 
such motions, y 1 and y 1 , which are termed objectively equivalent. 

From (3) 

y 1 = Q 1 j7 J + b 1 (7) 

Differentiating (7) with respect to time we find: 

v 1 = Q 1 , y j + Q 1 . 7 j + b 1 (8a) 

J J 

However the tensor transformation of 7 1 from y 1 to y 1 defines 
components v 1 given by 

T = (ayVay^) 7^ = Q 1 ^ 7^ (8b) 

i — i 

Thus velocity components referred to y and y are not related 

by a tensor transformation. We therefore conclude, not unexpectedly, 

that velocity is not an objective tensor. 

Necessary conditions for objectivity may be developed for 
tensor fields and functions of any order. Employing the two motions 
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y 1 and y 1 defined above, the following objectivity constraints are 
found for vectors q 1 , second order tensors t 1J and second order 
tensor functions f^ 1 J ) . Recalling that the transformation qV 
is a function of time we require the following to hold for all 
times. 

VuctohA : 

q 1 = Q 1 . q J (9) 

t 1J = t^ 1 (10) 

TtvUson Function!)*: 

f 1J (t k i) = Q 1 p Q ;i qf pq (Q r k Q s 1 t re ) (11) 

Strain and Strain Rate : The initial and deformed lengths of 

a differential material line segment are given by dS and ds 
respectively. 

dS 2 =. dX 1 dX J Gjj (12) 

ds 2 = dx 1 dx^ g. . (13) 

' J 

where g.. and G T1 are the metric tensors of the deformed and 
ij IJ 

undeformed coordinates. The change in length may be 


*Tensor functions satisfying constraints of the form (11) are termed 
hemi tropic functions of their arguments. 
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represented as 


ds 2 -dS 2 = dx i dxjg ij --dX I dX J G IJ (14) 

The differentials dx 1 and dX^ are related by the deformation 
mapping (1) 

dX 1 = X 1 .dx 1 (15) 

Therefore (14) may be written entirely in terms of the deformed 
coordinates x 1 

ds 2 -dS 2 = (g. j -X I ji X J J G IJ )dx 1 dx J ’ (16) 

The local deformation may therefore be described by the Almansi 
strain tensor e. . . 

* J 

5 0/2)<g tr x I , i x J ,jG IO ) (17) 

Similarly the local rate of deformation may be characterized by 
considering the time derivative of ds 2 in (13). 

d(ds 2 )/dt = 2 d-jjdx^dx J (18) 

d ij = (l/2)( Vi ;j + (19) 

The Euler deformation rate d. . provides a complete, objective 

J 

representation of the local non-rigid component of the velocity field. 
It is of particular consequence in the present analysis by virtue of 
its linear dependence on the velocity gradient. In the presence of 
local rotation the angular velocity of its principal axes is found 
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from the skew-symmetric velocity gradient, or spin tensor, u>. . 

* J 


<o k = e ki, i(jJ i j (20) 

“ij = 0/2)(v i . r v j;i ) (21) 

where e kl "J is the permutation operator. 

Almansi strain and the deformation rate are related as: 



d_ 

dt 


(x" 1 ,/ ,j e mn) 


(22) 


Stress and Stress Rate : The Cauchy stress tensor 0 1J referred 

J.L 

to the deformed configuration is the i tn component of traction 
on a surface normal to the coordinate direction. The symmetric, 
objective Cauchy stress provides a complete description of the loading 
state at a point in a deformed body. Components of traction t 1 on 
planes of arbitrary orientation having normal vector components 
are found as linear combinations of the a 1J 


t 1 = a ij ’vj (23) 

The time rate of traction t^ is found by differentiating (23), 
noting that the normal vector is itself time dependent. 

t 1 " = (a^-cr^v^ )v. (24) 

»P J 

The stress rate d. . in (24) is not objective. 

• J 

We must therefore seek an alternative objective characterization of 
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the time rate of stress for use in the constitutive equations of 
elasto-plastic flow. 

The objectivity constraint (10) is necessary but not sufficient 
to define a unique objective stress rate. Such tensors have been 
developed by a number of investigators including Jaumann [10], 

Truesdell [11], and Oldroyd [12]. As suggested by Prager [13] for 
use in analysis of elasto-plastic flow the present analysis employs 
the Jaumann rate. 

Consider a stress S 1J referred to the principal axes of the 
deformation rate (19). The rotation of these axes at a point in 
a deforming body is given by the spin u>. . of (21). The relative 

I J 

orientation of the rotating coordinates and a fixed frame is given 
by an orthogonal mapping Q 1 *., defined by 

J 

^j = Q 1 pQ p j (25) 

The stresses S 1J and o 1J are at all times related by 

sij = Q i pQ j q o pq (26) 

The time derivative of S 1J is objective. Differentiating (26) and 
transforming the result back to a fixed reference yields the 
symmetric objective Jaumann stress rate a 1J 

aiJ = Q 1 p Q j qS Pq (27) 

a"' J = a 1 ' J + a 1 " (28) 

mm 
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The J aumann AttieAA note ptiovUeA a meoAutie. ofi the tone, Hate. o£ AtfieAA 
a& Ae&n by an obAetiveti pafiticipating In flotation a de.&o tuning 

-i 

continuum at cvetiy point. 

A Jaumann rate may be constructed as above for any tensorial 
quantity. If the quantity is intrinsically objective the Jaumann 
rate and material derivative (6) are equivalent. 
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1.2 Constitutive Theory 


GeneAaZ ChaAa.cX.eA ofi the Then Ay 

A constitutive formulation is derived for elasto-plastic flow 
of metals undergoing finite deformation. The theory is derived as a 
generalization of an elasto-plastic flow theory appropriate for 
analysis of infinitesimal deformation (See Fung [14]). Its 
application is restricted to analysis of homogeneous isotropic metals 
undergoing quasi -static isothermal deformation. 

The character of the finite deformation constitutive formula- 
tion is dictated by the adoption of three fundamental characteris- 
tics of the infinitesimal theory. 

1. The flow mode, elastic or elasto-plastic, is dictated 
by the behavior of a scalar loading function f which 
is dependent upon current stress state and deformation 
history. 


2 . 


f = 0 ; f = 0 : the cases of loading, elasto- 

plastic flow, and neutral load- 
ing, elastic flow, may be 
identified. 

f < 0 : elastic flow or unloading 

The deformation rate is a pAioAt assumed to be separable 
into elastic and plastic components whose dependence upon 
stress is independently defined. 


(29) 
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3. The constitutive equations take the form of first order 
differential equations relating time rates of stress and 
deformation. 

The generalization to the finite case employs consistent 
interpretation of the time derivatives of stress and deformation 
appearing in the constitutive differential equations. The resulting 
formulation obeys the constraint of material objectivity. Deforma- 
tion magnitude is limited only by an assumption that recoverable, 
elastic, deformation is infinitesimal. 

In the following sections elastic and plastic flow modes are 
separately defined and then assembled to provide governing equations 
relating Jaumann stress rate and total deformation rate. 

Elastic Flow Mode : Following Green and Naghdi [15] an elastic 

strain component is defined as the difference between total and 
permanent deformation. Explicit representation of this recoverable 
strain component* in terms of problem kinematics is neither 

sought nor necessary. It is defined only as the time integral of 
d.jj( e ), the elastic component of the deformation rate tensor. 

Stress and elastic strain are assumed to be related through 
a positive definite strain energy density function W^ e ^. The 
energy density is defined such that 

a ij ’ = 3 w(e)/ 3e . (e) (30) 

* J 

*( e ) is not an index. 
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Since the strain energy is positive definite and the elastic strain 
is infinitesimal we may approximate as a quadratic function 
of ejj( e ). It may then be inferred* from (30) that the stress and 
elastic strain are related as 

0 ij = E «M eh (') (31) 

Taking a Jaumann time derivative of (31) and inverting the result, 
yields an objective relation between the elastic deformation rate 
and the Jaumann stress rate. 

d ij (fi) = M i jkl ° kl (32) 

For an isotropic material the constitutive compliance tensor 
has the simplest isotropic form for a fourth order tensor. 

M 1jkl - k [1/2 <9ik9jl + 911%) ~ T77 9 fJ 9 k] ] < 33 > 

The constants p arid v in (32) are the shear modulus and Poisson's 
ratio of classical linear elasticity. 

Plastic Flow Mode : Equations governing infinitesimal plastic 

flow of work hardening materials may be inferred from a hypothesis 

*This portion of the development follows that employed in Green's 
method for development of the generalized Hooke's law of classical 
linear elasticity. (See Eringen [8L chapter 5.) 
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first enunciated by Drucker [16]. It is postulated that plastic flow 
produced by application and removal of a self-equilibrated stress 
field is restricted by the rate* inequality. 

a 1J e..(P)>0 (34) 

* J 

icic 

Three characteristics of plastic flow theory are implied by (34). 

1. The loading function, f = const., is a convex closed 
figure in stress space. 

2. The plastic strain rate is normal to f = const. 

3. Plastic strain rate and stress rate are linearly related. 
Development of a similar formulation for finite deformation 

is impeded by the absence of a unique choice of an objective stress 
rate tensor for use in a generalized form of (34). In the present 
theory the Jaumann rate is chosen for its conceptual simplicity 
and obvious physical interpretation. The implications of Drucker 's 
hypothesis are postulated as operational characteristics of the 
finite theory. 

The loading function is taken to be of the form 

f = 4>(a ij ) - K(W(P)) (35) 

*For infinitesimal deformation these rates may be taken as partial 
time derivatives, for which objectivity is an unnecessary constraint 
since the deformed and undeformed states are indistinguishable. 

^Demonstration of these consequences of Drucker' s hypothesis is also 
provided by Naghdi (17). 
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In (35) <j> is the yield surface and K is a work hardening parameter 
determined solely by prior plastic work, W^'. 

Yield surface dependence upon current stress state is restricted 
to the invariants of the deviatoric stresses. Plastic flow inde- 
pendence of hydrostatic stress, and material isotropy are thereby 
guaranteed. This restriction also eliminates the Bauschinger 
effect since the yield surface will expand isotropically in stress 
space. Hence <f> in (35) is written 

<t> = <#>(^2*^3^ 

where J 9 = (l/2)s^sJ.- 

2 . J (36) 

J 3 • (i/3)s 1 j sa k s k 1 

The plastic component of the deformation rate is taken to 
be normal to the loading function. 

d^P) = A 9f/9a ij (37) 

The proportionality constant A in (37) may be found from a 
consistency condition for plastic flow given in (29) as 

f = 0 (38) 
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Taking the Jaumann time derivative of (35) obtains 


f = i-i < = (3+/3a 1 ^) o ij -K'W^ P ^ = 0 
K' = dK/dW^P) 

The Jaumann rate appears explicitly in (39) only for stress since 
material derivatives of all other quantities present are objective. 
The rate of plastic work |/ p ^ in (39) may be expressed as 

w(p) = a ij d .>) 

* J 

Substituting (41) and (37) into (39) one may solve for A in (37). 
The resulting expression for the plastic deformation rate is 

d (p) = (9j>/aa ij ) (3<fr/3o k1 ) ;kl 
^ [ K'(a*/aa rs )a rs I 

in which a linear relation between plastic deformation rate 

and Jaumann stress rate is apparent. 

Plastic flow is fully defined by (42). Practical use of this 

result requires a choice of an explicit form for the yield surface 

<f> and definition of the work hardening function K. To complete 

the formulation we employ a construction proposed by Swedlow [7 ] 

for the infinitesimal case. 

Dimensional consistency in (42) is established by assuming <j> 



to have the dimensions of stress. It is therefore plausible and 
convenient to consider it as an equivalent stress* r e q. We then 
may define an equivalent plastic strain rate d ^ conjugate 
to t in the sense that 

T eq d eq = = ° iJd ij^ ( 43 ) 

Substituting for d^P) in (43) leads after some manipulation to 
the equivalence 

K' - 0/-r eq )(i eq /d eq (P)) (44) 

This expression for K' provides means for its determination 
from experimental data. Defining a total equivalent plastic strain 
e e q(P) as the integral of d (44) may be written 

K' = 2 Weq <P)/T eq (45) 

where we have defined an equivalent plastic modulus y e q(P) 
dependent solely on the equivalent stress. 

* 0/2) d %q/ dc e5 P) (46) 

Derivation of the plastic modulus from uniaxial tensile test data 

is described in detail in Appendix I. 

*eq not indi ces 
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Introducing the equivalent stress and modulus in the flow equations 
(42) yields 



For convenient reference (47) is written as 

d ij (p) - N ijM skl < 48 > 

N^i is hemitropic function of the deviatoric stresses. The full 
flow equation is objective. 

Elasto-Plastic Flow : The total deformation rate is simply 

the sum of its elastic and plastic components. Assembling 
equations (32) and (48) 

d ij = d ij (e) + d ij (p) = E M ijkl + N ijkl ] £ kl (49) 

It is convenient to rewrite (49) as 

2pd 1j = B, jk , ;h 

B W S (,/2)[9 1^jl + 9n9jk]-^7 9 tj 9 k> < 50 > 

+ (p/p (p) )t (aT eq /3 ° 1J)(a W 3 ° l:1) 

^ ^ (9t /3o rs ) o rs 

L eq 

In (50) u is the shear modulus of linear elasticity. 
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y 


E/[2(l + v)] 


(51) 


The flow mode is controlled in (50) by the modulus ratio y/y eq ' p • 
For elastic loading and unloading the plastic modulus becomes 
infinite and (50) reduces to the elastic equations (32). 

Perfect plasticity is specifically excluded from the 
formulation. Hence 


y pn (p) t 0 (52) 

eq 

and (51) may be inverted. 

= P 1Jk ' d, , (53) 

kl 

In (53) 


pi jkl 5 


w[ ( 9 1k 9 jl +9 il 9 jk + ^ < 9 ij 9 ld )] 

kl 


- 2 yy 2 


(3t /ao’J )(3T eq /3a K1 ) 
1 +Y 2 (3T eq /3o mn )(3T e q / 3a rs )g nir g n? 


(54) 


Y 2 = T (y/y ^) / [ ( 3T ^ ri / 3a rs )a rS ] 
eq eq eq 

The tensor P 1 ^ is a hem'i tropic function of the deviatoric 
stresses. It possess the symmetries 

pi jkl _ pjikl _ pklji 


(55) 


(56) 
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and cannot be decomposed into elastic and plastic components. 

The inverse constitutive equations (53) are objective. The 
equations are expanded for planar flow in Apendix II. 

Examination of the constitutive equations reveals two 
critical features of finite elasto-plastic flow. 

1. The constitutive equations cannot be integrated to 
define relations between total stresses and strains 
except under very restricted conditions.* The stresses 
must be proportional and the deformation must be homogeneous. 

2. The total strain, found as the time integral** of d.., 
cannot in general be decomposed into elastic and plastic 
components. The elastic and plastic deformation rate 
components are defined with respect to the instantaneous 
configuration of a body. This configuration reflects the 
previous history of both elastic and plastic flow. Hence 
the total deformation reflects problem dependent coupling 
of elastic and plastic behavior. 

These facets of finite elasto-plastic deformation provide the 
motivation for the rate viewpoint adopted for the entire analysis. 

*The usual distinction between incremental and deformation theories 
of plasticity (Fung [14], p. 476]) is identified but must be 
extended to exclude local rotation. 

**See (22). 
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1.3 Equilibrium Theory 

The Concept ofi Rate Equtitb^Uum 

The rate nature of the constitutive equations for elasto- 
plastic materials dictates that their deformation be viewed as 
a time dependent flow process. Hence requirements of mechanical 
equilibrium must be applied not only to instantaneous states 
of a deforming body but also to the flow itself. 

Equilibrium equations governing total stresses and their 
rates of change are developed below. The derivation of the 
total stress equations provides a model for establishment of 
appropriate rate equations. 

The flow is taken as quasi -static thereby allowing inertial 
effects to be neglected. Body forces are also excluded. 

Total Stress Equilibrium ; The net load applied to a body 
B in static mechanical equilibrium must be zero. The net load 
is found as the integral of surface traction (23) over 3B the 
boundary of B. 
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Applying the divergence theorem* to (57) yields the volume 
integral over the body 



Since B is arbitrary, field equilibrium equations 

0 ^ . ^ = 0 in B 

s J 


may be inferred from (58). 

Equations (59) are the familiar stress equations of 
equilibrium in terms of Cauchy stress a 1J . These objective 
equations are valid irrespective of material constitutive 
behavior and deformation magnitude. 


♦Sufficient smoothness of o 1J in B is assumed. 



Stress Rate Equilibrium : Rate mechanical equilibrium of a 

deforming body requires that the time rate of net applied load 
be zero. The net load rate is given by the material derivative 
of (57). Existing stresses are assumed to satisfy (59). 

at / t1ds ’ / (i1 + < ’ 1JvP ;p'j ,ds = 0 
9B 3B 

t 1 in (60) is given by (24). The additional term derives from 
the time dependence of 3B. 

Applying the divergence theorem* to (60) yields 

f (cji J ' . . - aP - * • vJ.J dl/ = 0 

J 9 J sj 

B 

Again noting that B is arbitrary we infer field equations 

aij - a Pi . vi. = 0 in B 
;j »J ;p 

Satisfaction of (62) guarantees that given an equilibrated stress 
field equilibrium will be maintained in the presence of time 
varying loading. The equations are objective. 

♦Sufficient smoothness of a 1J , a , v. is assumed. 

l 



The stress rate equilibrium equations (62) have been derived 
by a number of previous investigators. Similarly to the total 
stress equations (59) they govern the stress field irrespective 
of deformation magnitude and material behavior. In the case of 
linear elasticity it may be argued (Hill [18]) that the first 
term in (62) dominates the second by the approximate order of 
magnitude ratio (E/o 1J ). Thus for elastic deformation of metals 
under moderate stress (ct 1 J«E) the rate equation may be approximated 
as 

a 1 ' J . • = 0 (63) 

9 J 

Since the total deformation is infinitesimal the time integrated 
equilibrium error introduced by dropping the additional term will 
be small . 

For elasto-plastic deformation, however, Rice [ig] notes that 
the dominance of the first term is diminished by the reduction in 
material stiffness, (y e qP/y)«l. The approximate form (63) may 
still be employed. for infinitesimal deformation but the equilibrium 
error will be larger, of more concern, and problem dependent. In 
the present analysis for finite elasto-plastic deformation the 
complete stress rate equilibrium. equations must be employed. 
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I . 4 Governing Equations of Finite Elasto-Plastic Flow 
The lyvitial-and Boundcmy-Value PtiobZm lot t Velocity 

The elements of a complete theory of finite elasto-plastic 
flow are now in hand. Assembly of the constitutive and equilibrium 
equations provides governing equations for the velocity field in 
a deforming body. Initial and boundary conditions admissible to 
these equations may be defined. The formulation reduces to previous- 
ly established results for infinitesimal deformation. 

In principle the governing equations may be integrated over 
space and time yielding solutions for complete stress and deforma- 
tion histories. In practice such integration is possible in closed 
form only in a limited number of simple cases. However, the quasi- 
Vi near nature of the velocity equations facilitates efficient 
numerical solution. The potential for obtaining such solutions for 
complex problems is a primary motivation for the present theoretical 
development. 

The Velocity Equilibrium Equations : The velocity field in a 

deforming body is chosen as the primary dependent variable thus 
guaranteeing flow compatibility. Equations governing this velocity 
field are found by assembling the constitutive and stress rate 
equilibrium equations. For convenient reference the equations are 
reiterated below. 
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Constitutive Equation 6 : 


s'J - P iJk ' d kl 


(64) 


Rate Eqdtttbatum Equation s: - crP 1 ' . • v J . =0 (65) 

5J >J 5P 

Writing the Jaumann stress rate in terms of the material rate 
a ij " the constitutive equations may be combined with the equilibrium 
equations to obtain: 


- oPi ;j v J 


= 0 


( 66 ) 


Using the definitions of the deformation rate (19) and spin (21) 

) 

tensors (66) may be written entirely in terms of velocity. 

(i/2)[g 1 n (g tk g 1 m -g t m g lk ) - g li g k m g t n -g ti g k n g 1 m ][a ,T,,I1 v t;1 ] ; . 

+ [R k 1 ij 'v j;1 ];i ♦° k Pv 1 '. pl - 0 


(67) 


This result, the velocity equilibrium equations, governs the 
instantaneous spatial dependence of the velocity field. Presuming 
knowledge of the stress field, these differential equations are 
linear at an instant of time. Furthermore they provide a quasi- 
1 inear model for the entire deformation process. Their solution 
involves simultaneous, but decoupled, integrations with respect to 
space and time. 


33 



Immediate integration of the velocity equilibrium equations with 
respect to time would provide equations governing the total deforma- 
tion. However such integration is possible only for homogeneous 
deformation under proportional loading, the same restrictions which 
limit such integration of the constitutive equations. Thus the rate 
viewpoint initially aidopted proves to be a viable approach not 
limited in application to specific classes of loading and deformation. 

Complete Problem Definition : Complete definition of the finite 

elasto-plastic deformation problem includes the governing equations 
(67), specification of material properties, and prescription of 
initial and boundary conditions for the dependent variables. 

HaXeAMil ?n opojvtiu,: The elastic constants defined by (33) must 

be known. Work hardening plastic flow character of the material 
is defined by the equivalent plastic modulus prescribed as a 
function of equivalent stress t 6C) . Any history of prior plastic 
deformation is reflected in this function. 

IviubLal CovuLUajoyi!>: The initial configuration. Bo, must be defined. 

Initial stress and velocity will normally be taken as null fields. 
However certain non-zero initial fields are admissible; the stress 
must be in equilibrium and the velocity single valued. 
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Bound&ny CondbtixonA : Admissible boundary conditions include 

prescription of both traction rate* and velocity on the boundary 
3B of the deforming solid. Traction rate and velocity vectors 
prescribed at the same point on 3B must be orthogonal. Three 
classes of problems may be identified. 

1. fundamental problems in which either traction rate 
or velocity are prescribed on the entirety of 3B, 

2. mixed problems in which traction rate and velocity 
are prescribed on distinct portions of 3B, and 

3. mixed-mixed problems in which traction rate and 
velocity are prescribed on the same portions of 3B 
and are limited by the orthogonality constraint cited 
above. 


♦Traction rate t is prescribed either explicitly, as in the case 
of pressure loading, or implicitly from knowledge of the total 
load rate, T 1 , applied to a finite portion, 3By, of 3B. In the 
explicit case from (24) 

t 1 ' = (o^-aPM.pjvj 

Implicitly the traction rate is given by T 1 through (60) 

V = J (t 1 + a ij vP. p Vj)ds 

3 By 
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The Limit of Infinitesimal Deformation : The velocity 


equilibrium equations (67) are valid irrespective of deformation 
magnitude. Reduction of these equations to forms previously 
established for infinitesimal deformation derives from the 
assumption that the stress-velocity gradient coupling terms 
make negligible contribution to the nature of the total deforma- 
tion. The distinction between deformed and undeformed coordinates 
becomes unnecessary. Under these assumptions (67) becomes 

[P kl %lLi “ 0 (68) 

These are the governing equations of infinitesimal elasto-plastic 
flow developed by Swedlow [7]. 

In the absence of plastic flow P^ 1J ' in (68) reduces to the 
constant linear elastic form in (31). The velocity 

equilibrium equations may then be immediately integrated with 
respect to time. The resulting equations (69) are the Navier 
equations for displacement u^ of classical linear elasticity. 

{[l/(T-2v)]g iJ g k p + g 1 ^} u 1;jk = 0 (69) 
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II. STRATEGY FOR PROBLEM SOLVING 


The solution to a problem of finite deformation must include 

• 1 / 

complete histories of the deformation mapping x n (X ,t) and the 
stress field a 1J (xk,t). Construction of this solution for the 
elasto-plastic case requires simultaneous integration of: 

1. the velocity equilibrium equations (67) to determine 
the velocity field in the time varying domain B, 

2. the velocity field with respect to time in order to 
determine B, and 

3. the constitutive equations (64) with respect to time, 
thus determining the stress field a 1 J(x* < ,t) in B. 

As has been previously noted analytic solutions may be found only 
for homogeneous deformation under proportional loading. Numerical 
solution is unavoidable for more general -problems. 

I 

Solution for the dependent variables as continuous functions 
of time requires an iterative approach. A variety of techniques 
might be employed including, for example, relaxation or predictor- 
corrector methods. While such iterative solution is feasible an 
enormous amount of computing might be anticipated with no guarantee 
of numerical stability. 

A more economic approach is suggested by the quasi -linearity 
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of the problem. Rather than seeking a continuous solution for the 
time varying configuration and stress field we restrict our attention 
to the behavior of these quantities at a finite number of times 
during the deformation. The total deformation is approximated as 
a sequence of incremental deformations. 

Adoption of the incremental viewpoint toward problem solving 
allows spatial and time integration to proceed sequentially rather 
than simultaneously. Spatial integration of the instantaneously 
linear velocity equilibrium equations provides the velocity field 
in B at time t. Subsequent integration of the velocity field and 
constitutive equations over a time increment 6t yields the con- 
figuration and stress field at a new time t + 6t. A new spatial 
problem for the velocity at t + St may then be defined. The 
computational efficiency of the incremental approach is immediately 
apparent. A complete problem is solved by a sequence of linear 
analyses. No iteration is required.* 

The availability of the complete theoretical rate formulation 
provides a distinct advantage for the incremental numerical 
solution procedure; Problem solving capability is not tied to 
particular numerical techniques. 

*A single exception exists. If elastic unloading occurs during an 
increment the plastic modulus becomes infinite and the analysis 
for that increment must be repeated. The iteration is, however, 
closed in the sense that it need be continued only until the local 
modulus value is consistent with the. behavior, of the point in 
question." ' ' " ’ ' ' " '■ 
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Irrespective of the numerical procedures chosen solution 
accuracy is controllable. The incremental model assumes that the 
velocity equilibrium equations written at an instant during the 
deformation provide an acceptable approximation over a finite 
time step. Note, however, that the coefficients in these 
equations are stress dependent and the configuration of the 
deforming solid defines their domain of integration. Hence the 
degree of approximation is dependent upon the variation in stress 
and configuration during a time step. The analyst retains 
control over the error in modeling a problem through his choice 
of time step size and is assured of convergence to a precise 
representation as the step size tends to zero. Thus increased 
time step size provides a less accurate solution at reduced 
expense and vice-versa, a measure of control not available 
in iterative solutions. 

The incremental approach to problem solving provides a 
vehicle for realization of the full potential of the rate 
formulation for finite deformation of elasto-plastic solids. 

No inherent restrictions exist upon loading type, geometry 
or deformation magnitude. The utility of the analysis is 
limited only by the availability of requisite material property 
data. 
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III. NUMERICAL SOLUTION PROCEDURE 


OvGAuXeu) 

The incremental approach is adopted for numerical solution 
of the finite elasto-plastic deformation problem. A finite 
element technique is employed to reduce spatial integration of 
the field problem of Section 1.4 to algebraic form. The overall 
procedure for analysis of a deformation increment is summarized 
below and developed in detail in subsequent sections. The 
procedures described are valid for analysis in three spatial 
dimensions in any coordinate system. 

The deforming solid is partitioned into an array of contiguous 
sub-regions or finite elements. Behavior of the complete solid 
is modeled by coupling these elements at a finite number of 
common points or nodes. Within each element spatial dependence 
of field -variables is approximated in terms of nodal values 
which become "the principal unknowns of the numerical problem. 

Linear algebraic equations governing the nodal velocities 
at the beginning of a time step are developed by applying the 
Galerkin method* [20] to the velocity equilibrium equations. 


*Note that for finite elasto-plastic flow energy principles are 
not available as a basis for finite element solutions. Hence 
an approach based solely on the governing differential equations 
is required. This point is further discussed in Section III.l. 
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These equations, termed the rate stiffness equations are of the 
form* 

T a = K«0V& (a, 6=1 M) (70) 

where Ta and V& are nodal loading rates and velocities, 
respectively. The range M is the total number of degrees of 
freedom associated with the finite element model . The rate 
stiffness in (70) depends upon the instantaneous configura- 
tion, prior plastic deformation, and the existing equilibrated 
stress field. Boundary conditions of the problem must prescribe 
precisely half of the 2M variables, load rates and velocities. 

Solution of the rate stiffness equations for the unknown 
nodal quantities provides the basis for evaluation of a deforma- 
tion increment. Nodal coordinates and loads as well as stress 
and strain fields in the elements are found by integration with 
respect to time. Values of these quantities at the beginning 
of the increment provide initial values for this integration. A 
new problem for the nodal velocities, may then be defined and 
the incremental procedure repeated. 


*Greek superscripts indicate matrix character. Repeated Greek 
superscripts are to be summed over an indicated range. 
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A single exception to the above procedure exists. Elastic 
unloading is signalled by a decrease in the equivalent stress 
in a region previously deforming plastically. When this occurs 
the rate stiffness must be recomputed to reflect elastic behavior 
in elements comprising that region at that time, and the incremental 
solution repeated. 

The incremental method described above has been implemented 
for analysis of planar deformation in plane stress or plane strain. 
Equations governing numerical solution of planar problems are 
given in Appendix III. 
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III.l Spatial Integration : The Galerkin Method 

The use of finite elements for spatial integration of the 
velocity equilibrium equations is contingent on the availability 
of a method for transforming these equations to an algebraic form 
involving a finite number of dependent variables. Discretization 
procedures developed for analysis of infinitesimal deformation, 
Zienkiewicz [21], provide a model for the present finite case. 

Reduction to a finite number of variables is accomplished 
by approximation of element field variables in terms of their 
nodal values. These nodal values are taken as the dependent 
variables of the numerical problem. Algebraic equations 
governing these variables must be derived from the governing 
differential equations. 

In analysis of infinitesimal deformation advantage is 
taken of the symmetric nature of the governing equations, e.g. 
the Navier equations of linear elasticity. Problems governed 
by such equations admit alternative statement as the variation 
of quadratic functionals. The variational problem may be 
extremal, e.g. minimum potential energy, or stationary, e.g. the 
Reissner [22] principle. The functional is written in terms of 
the finite element field variable approximation, integrated over 
each element and the variation taken with respect to the nodal 
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variables. The resulting linear algebraic equations are then 
solved for the nodal quantites. This entire procedure is known 
as the Ritz method. 

Application of the Ritz method to the present case of finite 
elasto-plastic deformation requires that the velocity equilibrium 
equations be symmetric. Writing these equations in terms of a 
linear differential operator L 1 , however, 

L 1 ’ {Vj} = 0 in B (71) 

and the boundary conditions as 

f 1 ( v j ) =0 on 3B (72) 

we find 


Jl L 1 '{Vj}]w i dl / 


J [L^Wj}] v. dl/ 


In (73) v. and w^ are independent, single-valued velocity fields 
satisfying (72). Symmetry requires that (73) be an equality. 
Thus the Ritz method is not applicable in the present case. An 
alternative integral method is required which, while admitting 
the use of a finite element field variable approximation, does 
not restrict the nature of the governing differential equations. 
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The present numerical solution procedure employs the Galerkin 
method which is applied directly to the governing differential 
equations and is appropriate irrespective of their character. Rate 
stiffness equations are developed for single elements and assembled 
to define master equations for a complete problem. 

The velocity field v^ in an element B is approximated* by v^ 

T k 

V. = V i = r«eVfy“(x K ) a,B=l N (74) 

where V e are nodal velocities, is dependent upon the nodal 

coordinates, is a vector of functions of x k and N is the 

number of degrees of freedom associated with the element. The 

(j>“ in (74) are prescribed functions of x providing an approximate 

representation of the spatial variation of v n - in B . The matrix 

1 m 

r ag is defined by requiring that evaluation of v. at the nodal 
positions yield V e . 

The Galerkin method** is based on the observation that if 
the in (74) are considered independent then requiring v^ to 
satisfy (71) as N tends to infinity implies orthogonality of 


♦Combined matrix tensor notation is employed. Greek superscripts 
denote matrix character and Latin indices denote tensor com- 
ponents. Thus the elements of <f>“ are first order tensor components. 
The overscript T denotes a matrix transpose. 

**This discussion is intended to communicate the essence of the 
technique. Rigorous discussions may be found in Kantorovich 
and Krylov [20] and Rektorys [23]. 
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( 75 ) 


each of <f>? (a=l,...,N) to L 1 {v n -} in B . Thus 
1 j m 

f [L 1 (Vj }.] *« dl/ = 0 a=l , . . . ,N 
B 

m 

The orthogonality conditions (75) provide N linear algebraic 
equations for the unknown of (74). Expansion of (75) for 
finite N yields the rate stiffness equations (76) for the 
element B m> The range of all Greek superscripts is N. 

y'Uaya t) r 6 V dl/ 

3B 

m 

_ 


(76) 



<f>. . 

1J 


(1/2) (+, , + 

1 »J 


j;i' 


(77) 


(78) 
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The surface integral in (76) corresponds to the time rate of load 
(See Section 1.4) on the element boundary and may be taken to 
define a vector of nodal load rates T a . Hence (76) is written 

jo = K ^V a n,a=l . ,N (79) 

The rate stiffness matrix K na is full; depends upon the element 
configuration through the nodal coordinates in r ag and upon the 
existing equilibrated stress state. In general the rate 
stiffness matrix is not symmetric. 

Having written (79) for each element, master rate stiffness 
equations are written for the entire body by summing the load 
rate components at each node. Thus follow the equations 

T a = K a6 V e a, 3=1 , . . . ,M ' (80) 

where the range M is the total number of degrees of freedom 
associated with the assemblage of finite elements used to model 
a complete body. The stiffness matrix K a & is M x M, sparse 
and not symmetric. Banded coefficient structure may be achieved 
in K 01 ^ by appropriate construction of the finite element map. 

Solution of the linear algebraic system (80) requires 
specification of M of the 2M variables, T a and V a . At internal 
nodes equilibrium requires T a to be zero. At boundary nodes 
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either T“ or V 01 may be prescribed* subject to the same ortho- 
gonality restriction governing the boundary conditions of the 
original field problem of Section 1.4. Solution of the equations 
may procede by any convenient technique and yields full knowledge 
of the velocity field in B and loading rate on 9B as character- 
ized by their nodal values. 


♦Boundary loading is prescribed either as rate of # total load T 
at nodal positions or through the traction rate t. on 3B. In 
the latter case the rate stiffness equations are Employed in 
the form (76) and the requisite surface integrations performed 
numerically. 
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III. 2 Time Integration 


The solution of the rate stiffness equations at time t 

O 

provides a basis for evaluation of ensuing changes in deformation 
and loading occurring during a small but finite time step <st. 

Time dependent variables are expanded in Taylor series about t , 
truncated to linear form and evaluated at t +<$t. The procedure 

O 

allows evaluation of the deformation, nodal coordinates and 
element strains, and loading, nodal forces and element stresses, 
at t +6t, the end of the time step. Requisite time derivatives 

o 

are provided by the rate stiffness equation solution at t . 

O 

The procedures of this section provide sufficient data, 
nodal coordinates and element stresses, for evaluation of the 
rate stiffness at t +6t, the beginning of the next time 

O 

step. Although more sophisticated time integration schemes 
might be employed, including higher order Taylor series 
representations, the linear technique proves adequate in practice. 

Nodal Coordinates.: The time dependent nodal coordinates 

are expanded as 

x a (t) = x a (t ) + [x a (t ) ] 6 t + ... (81) 
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Hence the nodal coordinates x a at the end of a deformation increment 
may be approximated at (t +<st) as 

O 

x a = x a + V a fit (82) 

O 

where x a are the nodal coordinates at t and V 01 are the nodal 

o 

velocities. Total nodal displacements U a are subsequently 
defined by 

U a = x“ - X 01 (83) 

where X° are the nodal coordinates of the original (t=0) 
undeformed configuration. 

Element Strains : The element deformation mapping function 

(1) at (t +<st) is written 

O 

xV. = X 1 + U 1 (84) 

where the total displacement field U* is taken as a function of 
the undeformed field coordinates X*. 

The element total displacement field is approximated in terms 
of nodal quantities in the same manner as was the velocity field 
in the development of the rate stiffness equations in Section III.l 

U I (X K ) = r ae U e <j> aI (X K ) (85) 
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Substituting (85) into (84) the element deformation gradient is 
written 

x 1 . = g 1 . + gV“ e u6 4> aJ > (86) 

il i O J 

Inverse deformation gradients axVax 1 are defined by observing 
that 

(axVsX 1 ) (axVaxJ) = g 1 . (87) 

J 

A variety of strain measures may be computed from the 
total deformation gradients (86). Almansi strain, for 
example, is given by 

e.. = (l/2)[g . - (sxVax 1 ) (3X J /3x j ) G.J (88) 

ij ij IJ 

It should be recalled that the strain is not separable into 
elastic and plastic portions. 

Nodal Loads : The time dependent nodal loads may be 

expanded as 

F a (t) = F a (t ) + [F a (t )] fit + ... (89) 

o o 

Truncating this series to linear form the nodal loads at the 

end of a time step are approximated as 

\ 

F a (t +fit) = F a (t ) + T a <St . (90) 

O O 

where T a are the nodal load rates of (80). 
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Element Stresses : Evaluation of the element stress field 

requires distinguishing between Cauchy stress cr 1J referred to 
the instantaneous configuration of an element and a Kirchhoff 
stress S 1 ^ referred to the configuration at t . Both o 1J and 

o 

S 1J are time dependent and are equal at t . Denoting the 

o 

o • * 

Cauchy stress at t as a J we define the Kirchhoff stress such 

O 

that 


lim 

t+t 


[S ij (x k ,t)] 



°k 

where x are the element field coordinates at t . At t +6t 

o o 

the two stress fields are related by 


(91) 


a 1 ^ (x k ,t + t) = (1/J) (8xVsx k ) (BX'Vsx 1 )S kl (x r ,t +6t) (92) 

o o 

where J = ) ax^ / ax j ] . Hence the final Cauchy stress may be 

found from the final Kirchhoff stress using the transformation (92). 

S' J is approximated at (t + t) by the truncated Taylor 
series 


S 1 3 (t +6t) = aiJ + S 1j 6t (93) 

O 

The Kirchhoff stress rate S 1J in (93) is computed from the 
Jaumann stress rate c 1 J at t . 


52 



(94) 


s'i 

where from (53) 


o'* + Pi d k k - 

K p 


n • • 

,iPdJ 


a ij = p 1Jkl d 


kl 


and the deformation rate d. . is computed from the element 

* J 

velocity field approximation (74). 

Computation of o 1J using (92) is facilitated by approximating 
the element incremental deformation mapping x (x J ) in the same 
manner as the total deformation of (86). Thus 

x 1 = x 1 + u 1 (x k ) (95) 

The incremental displacement field u 1 is represented in terms 
of incremental nodal displacements u e as 

T . 0 

u 1 ^) = r a ^u g ij) a1 (x k ) (96) 


where the incremental nodal displacements are determined from (82) 

u B = V e fit (97) 
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III. 3 The Limit of Infinitesimal Deformati on 


It has been noted in Section 1.4, that the velocity 
equilibrium equations reduce to well established results in the 
limit of infinitesimal deformation. It is therefore reasonable 
to expect similar limiting behavior of the algebraic equations 
governing finite element solutions. This expectation is reinforced 
by noting (Kantorovich and Krylov [20]) that application of the 
Galerkin method to symmetric differential equations, such as those 
governing the infinitesimal case, yields precisely the algebraic 
equations derived from a Ritz approximation. 

For infinitesimal deformation stress-velocity gradient 
coupling is assumed negligible and deformed and undeformed element 
configurations are taken to be indistinguishable. The rate stiff- 
ness equations (76) become 


T n = K nB v e 


where 


(98) 


Kna 



pijk! r 6ct } dV 


B 


m 


(99) 


The stiffness matrix K n “ is symmetric and positive definite and is 
recognized as that governing infinitesimal deformation, see, for 
example, Zienkiewicz [21] p. 16. 
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For elasto-plastic flow T a , V a in (98) are interpreted as 
increments of nodal load and displacement. Total strain (assumed 
«1) and stress are found as simple sums of incremental results 
since the reference frame transformations of Section I I I. 2 become 
identity operations. In the case of elastic behavior (98) 
integrates directly to provide linear relations between total 
nodal loads and displacements. 
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IV. EVALUATION OF SOLUTION CAPABILITY 


IV. 1 The FIPDEF Program 

The procedures of Section III have been implemented to 
provide incremental solution capability for problems of planar 
finite deformation. The computer program, FIPDEF (Finite 
Plastic DEFormation) , performs analysis of elastic and elasto- 
plastic materials deforming under conditions of plane stress 
or plane strain. Elastic unloading and subsequent reverse 
plastic loading are automatically treated. 

The program employs triangular finite elements defined 
by nodes at their vertices. The velocity field is approximated 
within each element by assuming linear spatial variation. Rate 
stiffness equations for these elements are developed in Appendix 
III. General program logic is shown in Figure (2). The FIPDEF 
program is written in Fortran IV and is operational on the 
Univac 1108 system at Carnegi e-Mellon University. 

Problem definition includes the initial finite element 
map geometry, material properties and incremental histories of 
nodal boundary conditions. Material behavior is described by 
the elastic constants, a proportional limit value of octahedral 
stress and pointwise specification of an octahedral stress- 
octahedral plastic strain curve.* 

*See Appendix I. 
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Figure 2 

GENERAL STRUCTURE OF THE FIPDEF PROGRAM 


1. Read and check input. 

2. Set initial elastic material properties. 

3. Generate rate stiffness matrix for the N th 
time step. (UNsNINC) 

4. Solve for nodal velocities and load rates for 
step N. 

5. Evaluate element stresses at the end of step N. 

6. Check all elements for load reversal. 

a. if reversal has occurred in any element(s); 
modify element material properties accord- 
ingly and return to 3. to repeat the N th step. 

b. if loading continues in all elements evaluate 
their material properties for step N+l . 

7. Evaluate nodal loads and coordinates and element 
strains at the end of step N. 

8. Output 

9. If N < NINC: go to 3 for step N+l. 

10. If N = NINC: STOP. 
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IV. 2 Verification Analysis 


Evaluation of the finite deformation analysis capability 
provided by the FIPDEF program is difficult since comparison 
with analytic solutions is possible for a limited number of 
problems. In particular analytic solutions may be obtained only 
for homogeneous deformation under proportional loading. Despite 
their conceptual simplicity, however, such problems retain 
considerable nonlinearity and therefore provide an acceptable, 
albeit limited, basis for assessment of numerical results. 

Furthermore, since the deformation of individual finite 
elements employed in FIPDEF analysis is restricted to homogeneous 
form, the behavior of these elements may be completely evaluated. 
Verification of the accuracy of individual element response is 
significant for two reasons. 

1. It demonstrates the viability of both the rate theory 
and the incremental approach for solution of finite 
deformation problems. 

2. Accurate prediction of individual element behavior is 
a necessary condition for accurate solution of inhomo- 
geneous deformation problems requiring the use of 
element arrays. 
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It should be recognized that the use of element arrays to 
solve inhomogeneous problems may of itself introduce error in 
numerical solutions. Quantitative assessment of this discre- 
tization error is not possible on the basis of the homogeneous 
deformation problems for which analytic solutions are available. 
Some indication of the potential significance of this type of 
error can be inferred from the numerical solutions of planar 
necking problems which are presented and evaluated in Section V. 

FIPDEF and analytic results are compared below for three 
classes of homogeneous finite deformation problems: extension, 

simple shear and simultaneous extension and rigid rotation. 
These three types of problems span the range of possible finite 
deformations. The non-rigid portion of any continuum motion 
may be viewed as the simultaneous occurrence of dimensional 
changes, such as those associated with finite extension, and 
combined shearing and rotation of the form found in simple 
shear. The third problem, extension and rotation, provides 
a means for verifying the objective character of the analysis 
since the solution must be tensorially independent of rotation. 

Complete solutions for the verification problems are 
developed in Appendix IV. Salient features of these solutions 
are discussed below and compared with the results of numerical 
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analysis. It is demonstrated that accurate, objective numerical 
solutions can be obtained for deformations involving dimensional 
changes of a full order of magnitude and rotations of up to 45 
degrees. These upper limits do not reflect deterioration of 
the numerical results but rather the judgement that most 
deformations of practical interest will be within these bounds. 
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IV. 3 Finite Extension 


The unit cube of Figure (3a) is deformed into a bar, Figure 
(3b), by prescribing its horizontal or x dimension as a function of 
time. The deformation at any time is fully described by the 
coordinate stretch ratios 


V E Ijl 

A X o 


x y " l y /l o 


( 100 ) 


= in 

Z Z o 

in which l 0 is the original unit dimension. No shearing deformation occurs, 
The deformation is produced by uniaxial loading (a ^0, 

A 

ay=0) and may procede under conditions of plane stress or plane 
strain. In plane stress we require 


while in plane strain 


°z = 0 


A z - 1 


In either case the applied load is found as 


P v = a A 
x xx 


( 101 ) 
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where 


°x = °x (e x) 


e x ^n X x 


and the cross-sectional area A is found from 

A 


A = A (x ) = X X 
x x' x y z 


In plane strain we have the additional result 


a z = a z (x x> 


( 102 ) 


(103) 


(104) 


Full solutions for load, stresses and area are developed in 
Appendix IV for plane stress and plane strain extension of elastic 
and elasto-plastic materials. 

For purposes of verification analysis elastic materials are 
defined as those governed by (32) irrespective of deformation 
magnitude. To facilitate analytic solution elasto-plastic problems 
employ the bilinear material property representation of Figure (4).* 
While the validity of material property models restricts 
engineering application of the FIPDEF program, it does not 
affect evaluation of solution accuracy. The only concern here 
is consistent definition of verification problems to be solved 
numerically and analytically. 


*The effective stress and plastic strain of Figure (4) are directly 
related to the "equivalent" quantities of Section 1.2 (see Appendix 
I). 
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FIPDEF and analytic solutions have been compared for problems 
of plane stress and plane strain extension for elastic and elasto- 
plastic materials. All FIPDEF analyses employ the two element 
model of Figure (3c). The deformation was prescribed in terms of 
incremental x displacements of the appropriate nodes in this 
model. With one exception all displacement increments were one 
percent of the original unit dimension of the cube. For elasto- 
plastic problems plastic flow was first established by 
prescription of a number of very small ( « .01 ) incremental 
displacements. The increment sizes were arbitrarily chosen. 
Detailed study of increment size effects was not undertaken. 

Agreement between analytic and numerical results was 
excellent in all cases. Applied load, stress and deformed 
cross-sectional area were predicted within one percent of 
the analytic results. 

Figure (5) is a plot of P , a , and A as functions of X 

AAA X 

for an elastic material in plane stress. It illustrates a 

distinguishing characteristic of all problems of finite extension. 

The applied load-stretch relation is multiple valued;* a maximum 

load P is attained at a critical, property dependent, stretch X . 
c c 

*The effect is observable in the data from any tensile test. 
Discussion of its physical significance is undertaken in Section 
V. 
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The maximum load and critical stretch may be analytically predicted 
by solving the equation 


d ( p X > /d(x X ) =( d V dX x )A x + a x (dA X /dx x ) = 0 

The plausibility of the maximum load phenomena is evident in 
(105) since, for the materials considered, (dA /dx ) < 0. 

X X 

Analytic expressions for P c and X £ are given in Table ( I V - 1 ) , 
Appendix IV. 

Numerical comparison of FIPDEF predictions and analytic 

results are given in Tables (1) and (2) for elastic and elasto- 

plastic problems respectively. Critical stretch, maximum load 

and critical axial stress are predicted within one percent. All 

results are dependent upon elastic and plastic properties. 

Sensitivity to the elastic Poisson ration is shown in Tables 

(1) and (2) to illustrate the resolution which is attainable. 

Load-stretch results for plane stress and plane strain 

extension of elastic and elasto-plastic materials are given in 

Figures (6) and (7) respectively. The results are normalized on 

the critical values P„ and X . FIPDEF results and analytic 

c c 

predictions coincide over the full range of deformation considered 


(105) 
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Table 1 


Homogeneous Extension of Elastic* Bodies 
Comparison of FIPDEF and Analytic Results 
at Maximum Load 


Plane Stress: 


n * l ^ U V* 

Poisson s _c 9 f, 

Ratio (lbs x 10 ) (lb/IrT x 10" D ) 

FIPDEF Analysis FIPDEF Analysis FIPDEF Analysis 


0.3 

6.15 6.13 

5.27 

5.30 

1.67 

1 .67 

0.4 

4.62 4.60 

3.48 

3.49 

1.25 

1.25 

0.5 

3.70 3.68 

2.71 

2.72 

1.00 

1.00 

Plane Strain: 






0.3 

9.46 9.43 

10.16 

10.29 

2.57 

2.56 

0.4 

6.59 6.57 

4.46 

4.48 

1.79 

1 .79 

0.5** 

4.93 4.91 

2.71 

2.72 

1.34 

1.34 

*Note: E = 

10 6 lb/in^ in all 

cases. 





**Poisson's ratio of 0.5 is inadmissible in plane strain analysis 
(See Appendix II). Since it corresponds to elastic incompressir 
bility the present analyses were performed using a high ratio 
of bulk modulus to Young's modulus k/E = 10^, which approximates 
.'such material . : ~~ .... 
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Table 2 


Homogeneous Extension of Elasto-Plastic* Bodies 
Comparison of FIPDEF and Analytic Results 
at Maximum Load 


Plane Stress: 


Poisson's 

Ratio 

(lbs 

V 

x 10‘ 4 ) 

X c 

°c 

(lb/in^ x 

10' 3 ) 

0.3 

7.78 

7.72 

1.27 

1.27 

94.74 

94.20 

0.4 

7.65 

7.58 

1.25 

1.24 

93.29 

92.56 

0.5 

7.52 

7.44 

1,23 

1.22 

91.45 

90.91 


Plane Strain: 






0.3 9.29 

9.17 

- 1.44 

1.43 

127.6 

127.1 

0.4 9.20 

9.09 

1.42 

1.41 

125.1 

124.0 

0.5** 9.03 

8.92 

1.36 

1.36 

121.7 

121.0 


*Note: E = 10^ lb/in^ _ 

3 = da/devP) = 10^ lb/in for a > a, 
a = 8 x 10 4 lb/in 2 


**Poisson’s ratio of 0.5 is inadmissible in plane strain analysis 
(See Appendix II). Since it corresponds to elastic incompressi- 
bility the present analyses were performed using a high ratio 
of bulk modulus to Young's modulus k/E = 10 2 , which approximates 
such material. 
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IV. 4 Simple Shear 


Consider a uniform elastic continuum deforming according 
to the prescribed velocity field 

v x = 2ky 

v = v = 0 

y z 

The resulting deformation is homogeneous simple shear in x-y 
planes of the material. Thus a unit cube of material, Figure 
(8a), becomes an oblique prism of unit depth. Figure (8b). The 
deformation is completely described by the shear angle y in 
Figure (8b). 

y = tan - "* x 
x = 2kt 

The solution for the time dependent stress field resulting from 
this deformation is developed in Appendix IV. The non-zero 
stress components are 


a = v sin x 
xy 

a x = y (1-COS x) 
a = y (cos T-l) 


(106) 


( 107 ) 


( 108 ) 
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FIPDEF analysis was performed in both plane stress and plane 
strain employing the two element model of Figure (8c). Simple 
shearing deformation was developed by prescribing incremental x 
displacements of the upper nodes of the element map. All other 
displacements are null. For the unit dimension model these 
incremental displacements correspond to increments of x in (107). 
The increment size was 0.01. 

Analytic and FIPDEF stress results are compared in Figure 
(9). Shear stress is predicted within one percent and normal 
stress within five percent over the deformation range 0 s x s 0.7. 
Results from plane stress and plane strain analyses are identical. 
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IV. 5 Simultaneous Extension and Rotation 


Consider the problem of simultaneous rotation and unidirectional 
extension of a bar of elastic material. As shown in Figures 
(10 a,b), the x‘ dimension of the bar is prescribed as a function 
of time while its y' and z' dimensions are maintained at their 
original unit values. The bar simultaneous undergoes a rigid 
rotation e(t). This motion may be described in two objectively 
equivalent ways. In the primed, rotating coordinate frame 
only the extension process is observed. In the unprimed, fixed 
frame both extension and rotation must be taken into account. 

Components of stress and surface traction corresponding 

to the above homogeneous deformation are readily found by 

analyzing the problem in the rotating reference frame of Figure 

(10). This solution is developed in Appendix IV. It predicts 

non-zero normal stresses (a v ■, a wl , a, ( ) and normal tractions 

a y z 

(t x i , tyi) on x' and y 1 coordinate planes respectively. 

FIPDEF analysis, however, proceeds in a fixed non-rotating 
reference frame. Solutions developed in such a reference frame 
must be independent of rigid rotation, i.e., the analysis must 
be objective. Hence the FIPDEF solution must be related to 
the rotating frame solution by the objectivity conditions of 
Section 1.1. Denoting the results of fixed frame analysis by 
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unprimed quantities we require that the stresses be related as 


a = cr t cos 2 e + a i sin 2 e 
x * y 

a = a , sin 2 e + o , cos 2 e (109) 

y * y 

a z = a z . 

°xy = a y' ^ cose sine 

and the in- plane nodal forces as 


F ^ = F , W cos e - F sin e 
x -* y' 

F W = F. M sin 0 + F , (N) cos 0 

y x y' 


(no) 


In 009,110) e is the time dependent orientation of the extending 
bar. In (110) F M is the x component of total load at the 

A 

node of a finite element model. F is analogously defined 
in the x' system but corresponds to the surface tractions of the 
rotating frame solution. 

The elastic extension-rotation problem was analyzed 
using the two element model of Figure (10c). Incremental displace- 
ments were prescribed at all nodes such that the body underwent 
simultaneous increments of extension and rigid rotation about 
node 1 in that figure. Extension increments were one percent of 
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the original unit length. Rotation increments were 0.005 radians. 
The final stretch ratio, A x i , was 2.6 and the final orientation 
angle was 0.8 radians. 

Results obtained agreed with equations (109) and (110) 
within one percent over the full range of deformation considered. 
FIPDEF stress results are plotted in Figure (11) as functions of 
X i and e. Representative nodal force results are given in Figure 

A 

( 12 ). 

On the basis of these results, as well as those of similar 
analyses for plane stress, it is concluded that affects of finite 
rotation are properly treated in the FIPDEF program. The ob- 
jective character of the underlying rate theory is preserved. 
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Figure 11 


Stress Components: 

Combined Extension and Rotation 



1.4 

1.8 

Stretch X x , 

2.2 

2.6 

0.2 

0.4 

Rotation 6 (radians) 

0.6 

0.8 
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V. NECKING IN FLAT TENSILE BARS 


The numerical solution capability which has been developed 
for problems of two dimensional finite deformation of elasto-plasti 
materials provides a basis for the investigation of necking in flat 
tensile bars. The physical aspects of the problem as well as a 
number of significant previous analyses are discussed below. 
Subsequently the numerical problem is described and results of 
analyses under conditions of plane stress and plane strain are 
presented and discussed. The results provide some insight into 
the mechanics of the necking process as well as a basis for 
quantitative evaluation of the utility of stress-strain relations 
inferred from post-instability tensile test data. 



V.l The Physical Problem 

Consider a flat metal bar. Figure (13), initially of length 
l , width w and thickness t . In a tensile test the bar is 

o o o 

quasi -statically extended in the y direction by prescribing its 
deformed length (l > l ) as a function of time while the lateral 

o 

surfaces of the bar remain traction free. The present discussion 
is limited to tensile testing of metals at room temperature 
under atmospheric pressure. • 

Under the above conditions the bar will initially undergo 
a process of homogeneous extension requiring a monotonically 
increasing applied load.* For extensions in excess of some 
critical value the applied load is observed to decrease with 
increasing extension and subsequently a highly strained local 
necked region develops in the bar. Extensions corresponding 
to maximum load and necking initiation are found to vary with 
the material and to a lesser degree with initial bar geometry. 

As overall bar length continues to increase plastic straining 
continues in the necked region while previously yielded material 
in the remainder of the bar unloads elasticaVly. For convenience 
we define the onset of necking as that point at which elastic 
unloading first occurs. 

*The upper and lower yield point phenomenon characteristic of 
mild steel is excluded from this discussion. 
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Neck development in flat bars is observed to be dramatically 
dependent upon initial bar geometry. Nadai [2] reports that for 
thin bars (w /t > 10) necking takes the form of local thinning 

o o 

in oblique bands parallel to the x-y plane of Figure (13). For 
thicker bars (w /t < 7), however, necking appears as a symmetric 

o o 

reduction in width in an x-z plane of that figure. In bars of 
intermediate aspect ratio (7 < w /t < 10) complex combinations 

o o 

of these necking modes are observed. 

Experimental investigation of the necking instability 
phenomenon has not provided a sufficient basis for identification 
of criteria for either necking initiation or location of the 
necking region in an initially prismatic bar. Criteria for 
attainment of maximum load*, such as those employed in verifying 
the FIPDEF program (Appendix IV), provide only a lower bound 
for necking initiation. Experimental and very limited analytical 
results indicate, however, that necking occurs somewhat later 
than maximum load. The additional extension required is dependent 
upon the properties of the material tested and possibly upon 
initial bar geometry. 

Presuming the existence of a stress- and/or strain-based 

*The earliest prediction of maximum load criteria was provided 
by Considere [24] who considered simple extension of incompres- 
sible materials. 
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necking initiation criterion neck location might be determined 
by the distribution of small imperfections in geometry or 
material properties in seemingly perfect bars. Satisfaction of 
the necking criteria in the vicinity of one such flaw, or in a 
region of high flaw density, will cause additional prescribed 
bar extension to concentrate in the flaw region, thereby 
developing a neck. It has been shown that the intentional 
introduction of geometric flaws can be used to control neck 
location. Standard test procedures (ASTM [25] standard E-8) 
for flat bar specimens allow a local reduction of area of up 
to 10 percent for this purpose. It must be noted, however, 
that the above interpretation of necking initiation has not 
to date yielded complete understanding of the mechanics of 
the problem. 

The present investigation is limited to consideration of 
symmetric necking in flat bars of elasto-plasti'c material. In 
lieu of analysis in three spatial dimensions, the present effort 
is restricted to consideration of the limiting cases of plane 
stress and plane strain deformation in the x-y plane of Figure 
(13). Neck location is controlled by considering an initial 
geometric flaw in the form of an 0.5 percent area reduction in 
the plane y = 0, which becomes the root plane of the neck. The 



The analysis is further simplified by assuming the transverse 
tractions on the ends of the bar (y = ±1/2) to be zero. 
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V.2 Previous Analyses 


Tensile bar necking has been analyzed by a variety of 
techniques. The utility of these analyses has been limited 
largely by their consideration of idealized models of material 
behavior as well as, in some cases, an inability to predict both 
deformed geometry and internal stress and strain fields throughout 
the necking bar. A number of these analyses are described below 
to indicate the range of previous efforts as well as to provide 
a basis for comparison with the results of the present analysis. 

Bridgeman [1] developed a solution for the stress distributions 
in the root planes of necks in round and flat tensile bars. The 
analyses consider a material obeying a Mises yield criterion 
and neglect elastic deformation. An assumption of uniform 
strain in the root plane is made in both cases. The assumption 
is corroborated for the axisymmetric case by the experimental 
work of Bridgeman himself as well as that of Davidenkov and 
Spiridonova [26], No experimental evidence exists to support 
the uniform strain assumption in the planar case. The solutions, 
while physically plausible, are of limited use since a pfviohA. 
knowledge of neck geometry is required and no relationship is 
established between overall bar extension and development of 
the neck. 
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Full field solutions for tensile necking in flat bars have 
been obtained by a number of investigators through analysis of 
plane strain extension for rigid-plastic materials. The 
perfectly plastic case has been considered by Richmond [27] and 
Onat and Prager [28] and hardening plastic by Cowper and Onat [6]. 

These solutions may also be distinguished by the manner in 
which necking deformation is introduced. Richmond considers the 
prismatic bar as the limit of a V-notched bar and thereby has 
in effect introduced an initial imperfection. Onat and Prager 
explicitly consider initial shallow longitudinal grooves. Cowper 
and Onat, on the other hand, develop admissible solutions correspond- 
ing to incipient necking by demonstrating bifurcation of the 
solution to the homogeneous extension problem for a hardening 
material. 

Neck geometries are predicted by all three analyses. 

Richmond predicts a boundary profile which is convex at the root 
plane and tends to a linear profile away from the neck region. 
Qualitative agreement with limited experimental data is 
demonstrated. The Richmond profile, while differing drastically 
with the V shaped neck predicted by Onat and Prager, is qualita- 
tively similar to that shown to be admissible for hardening 
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material by Cowper and Onat. This similarity perhaps indicates 
that analysis of necking employing slight initial geometric 
imperfections provides results indicative of behavior of perfect 
bars. 

Analysis of necking has recently been extended to include 
consideration of elasto-plastic (work-hardening) material 
behavior by Chen [29] and Needleman £30]. Both investigators 
have developed numerical solutions for the axi symmetric case. 

The two studies are significant in the context of the present 
effort not only by virtue of their consideration of elastic 
deformation but also because both analyses indicate that necking 
initiation is distinct from attainment of maximum load. The 
observation is of increased significance since Chen considers 
an initially imperfect bar while Needleman establishes necking 
deformation using a bifurcation technique. 
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V. 3 The Numerical Problem 

Two numerical analyses have been performed, one in plane 
stress, the other in plane strain. Initial geometry, material 
properties and in plane boundary conditions were identical in 
both analyses. 

The undeformed bar of Figure (13) is represented in two 
dimensions by the finite element model of Figure (14). The 
model employs 600 finite elements defined by the positions of 
341 nodes. The finite element model is bounded by symmetry lines 
at x = 0 and y = 0 and represents one quarter of a complete bar. 
An initial geometric imperfection is introduced by reducing the 
cross section at y = 0 by 0.5 percent. The amplitude of the 
imperfection is shown greatly magnified in Figure (14). 

The bar is taken to be of isotropic homogeneous material 
whose inelastic deformation is described by the constitutive 
formulation of Section 1.2. Elastic deformation is characterized 
by the constants 

E = 10 x 10 6 lb/ in 2 
v = 0.3 

Work hardening plastic deformation is controlled by the effective 
stress-effective plastic strain relation plotted on both linear 
and logarithmic scales in Figure (15). Initial yield occurs for 

a e f = ay = 35 x 10^ lb/in 2 
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y 


Figure 14 Flat Tensile Bar: Finite Element Model 
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Thereafter the stress-strain relation is 


e (P) 

e ef 


A[(o f/ay)-!]' 


(a ef - K) / B 


0 < e^< 0.05 


e^> 0.05 
ef 


(in) 


where 


A = 5.0814 
N = 3.0 


K = 40 x 10 3 lb/in 2 
g = 50 x 10 3 lb/in 2 


The effective modulus ratio E /(da^/deg^ 1 "^ ) therefore varies 
smoothly from zero at the yield point to a constant value of 
200 for effective plastic strains in excess of 0.05. The 
plastic stress-strain curve is input to the analysis in the form 
of data points corresponding to (111). The spacing of these 
points in strain varies from 6e^= 10" 5 at initial yield to 
<$e^P)= 0.05 in the linear portion of the curve. 

The initial stress-free configuration of the two dimensional 
quarter symmetry model of the bar is shown in Figure (14). The 
initial length to width ratio is 3:1. The model is deformed by 
prescribing the history of incremental boundary conditions. Within 
each increment these boundary conditions are as follows: 

1. Normal displacements and tangential traction are zero 
on the symmetry boundaries x = 0 and y = 0. 
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2. The lateral boundary initially at x = w /2 is traction 

o 

free throughout the analysis. 

3. Extension is imposed by prescription of uniform positive 
normal displacement of the boundary initially at y = l /2. 

o 

Tangential traction on the boundary is zero. 

The magnitude of the incremental extensions was varied over the 
extension history. Initially small steps {6l/l z 10~ 3 ) were taken 

o 

to establish plastic flow over the entire bar. Incremental 
extensions were then gradually increased according to the algorithm 
6l/(l-l 0 ), = 0.05 until necking initiated. Subsequent incremental 
extensions were maintained at approximately 0.5 percent of the 
length of that portion of the bar which continued to deform 
plastically. The final overall stretch ratio l/l considered 

o 

was 1.43 in plane stress and 1.62 in plane strain. 

Approximately 130 loading increments were employed in each 
case, each increment requiring an average of one minute of 
computing time on a Uni vac 1108. The analyses were terminated 
when deformations were developed which were judged sufficient 
to permit reasonable assessment of the necking process. The 
results do not suggest any breakdown in the analysis. There is no 
reason to suspect that the analysis could not have been extended 
indefinitely although in the absence of quantitative criteria for 
prediction of tensile fracture such effort was not warranted. 
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The analyses provide full histories of deformation and stress 
associated with the process of extension and necking in plane 
stress and plane strain. The deformation is given explicitly by 
the deformed configuration of the finite element model. The 
stress field is represented in terms of its component values at 
the centroid of each finite element. 
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V.4 Necking in Plane Stress and Plane Strain 

Overall tensile bar behavior may be characterized in terms of 
relationships between extension and applied load provided by the 
numerical analysis. These relationships are discussed below and 
compared with the results of homogeneous analysis in which necking 
is not considered. Inspection of neck geometry and internal 
stress field histories provides some insight into differences 
between load-extension relations predicted for plane stress and 
plane strain extension. 

Load-Extension Response : Numerically established relationships 

between applied load and extension or engineering strain e" = &1/1 , 

o 

are given in Figure (16) for both plane stress and plane strain. 
Results are shown from the inhomogeneous analysis of Section V.3 
as well as from homogeneous analysis in which necking is disallowed. 
The inhomogeneous and homogeneous numerical analyses are distinguished 
only by the use in the latter case of a two member finite element 
model of an initially prismatic bar. The homogeneous approach 
provides what might be termed fundamental solutions for finite 
elasto-plastic extension which provide a convenient reference in 
discussion of the necking process. 

A maximum load phenomenon is evident in all cases shown in 
Figure (16). The critical extension e Q at which maximum 
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Applied Load P (lbs x 10 


Figure 16 


Load-Extension Data; 

Plane Stress and Plane Strain 




load occurs is significantly larger in plane strain than in 
plane stress. In both cases ^ for the perfect bar is slightly 
larger than that for the initially imperfect bar. 

These differences in critical extension may be explained 
qualitatively and in an approximate sense quantitatively on the 
bases of the Considere [24] criterion for attainment of maximum 
load. By solving the equation 

dP y /d^ =0 (112) 

for an incompressible material it may be inferred that at 
maximum load 

a y = e y = VS 013) 

The overscript bar denotes deformation measures averaged over the 
entire bar as opposed to local values * In (112, 113) P is the 

y 

tensile load, X = 1 + e = l/l is the uniform stretch of the 
y y o 

bar, o is the so-called true stress (load divided by current 

y 

area) and 7 = InX is the logarithmic or natural strain. The 

y y 

assumption of incompressibility, which is not made in the numerical 
analyses, is later shown to be a reasonable approximation for 
some purposes. Noting the linear nature of the stress-strain 
curve for large plastic strains. Figure (15), it may be predicted 
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using (113) that the maximum load natural strain in plane strain 
is approximately 1.7 times the corresponding strain in plane stress. 
The data of Figure (16) show the numerically established ratios 
to be 1.57 and 1.6 for the perfect and imperfect bars, respectively. 
The deviation between the above prediction and the numerical 
results is a consequence of the elastic dilatation considered in 
the numerical analysis. 

The difference in critical extension between the perfect 
and imperfect bars is qualitatively explained by the presence 
of the initial imperfection. In the imperfect bar the stress 
at the location of the initial area reduction will be somewhat 
higher than that existing in the perfect bar for the same overall 
extension. Hence (113) is satisfied for a slightly smaller (-10%) 
extension of the initially flawed bar. The difference in initial 
.minimum cross-sectional area similarly explains the development of 
lower loads in the imperfect bar than in the initially prismatic 
bar prior to attainment of maximum load. 

Post-maximum load behavior in both plane stress and plane 
strain is characterized by a more rapid reduction. in applied load 
than that demonstrated by the fundamental solutions of Figure (16). 
The load reduction in plane stress is noticeably greater than the 
reduction in plane strain. Six stages in the load reduction process 
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are noted in Figure (16) for convenient reference in discussing the 
associated necking deformation. The first stage is the maximum 
load point. 

The difference between plane stress and plane strain load 
reduction is more readily apparent in Figure (17a) in which post- 
maximum load response is plotted normalized to the critical 
values, maximum load P and critical overall stretch X . The 

C V 

fundamental solution is also shown for an incompressible 
material. This solution, identical for plane stress and plane 
strain, is found as 

P/P c = (T/x) [ln(X/X c ) + 1] (114) 

A more rapid deterioration of applied load in the case of 
plane stress extension is apparent in Figure (17a). This 
observation is consistent with the more rapid concentration of 
stretch at the root- plane of the neck in plane stress as shown 
in Figure (17b). The figure shows X, the average root plane 
stretch, as a function of X, the stretch of the entire bar. 

The data are again normalized to their values at the maximum 
load point, stage 1. The root plane stretch X is computed from 
the numerical results as (assuming incompressibility) 

X = A /A (115) 

O 
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where A (A) is the original (current) root plane area. 

O 

The data of Figures (16, 17) also suggest that necking does 
not initiate until somewhat after maximum load is attained. This 
observation has also been made by Chen [29] and Needleman [30] 
based on analyses of the elasto-plastic axisymmetric case. The 
present results also indicate that necking initiation, signified 
by the occurrence of elastic unloading away from the root plane, 
occurs somewhat later in plane strain than in plane stress. The 
normalized overall stretch ratio T/T c at necking initiation is 
approximately 1.03 in plane stress and 1.06 in plane strain for 
the material and initial geometry considered in the present 
analyses. 

Necking Deformation : The physical character of the necking 

process is illustrated by the deformation histories depicted in 
Figure (18). Configurations of the finite element model are shown 
which correspond to the undeformed state and the six stages in 
the necking process previously identified in Figure (16). Results 
are shown for both plane stress and plane strain. 

The geometric imperfection present in the undeformed state 
appears slightly amplified in the maximum load configuration, 
stage 1. Shortly thereafter necking initiates and subsequent 
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prescribed bar extension is increasingly concentrated in the vicinity 
of the initial area reduction, stages 2-6. Simultaneously the 
exi stance of decreasing load in combination with little or no 
reduction in area produces elastic unloading in regions removed 
from the neck. Unloading first appears in the center (x = 0) of 
the bar at the furthest boundary (y = U 2) from the initial flaw 
and thereafter spreads down the bar. The elastic-plastic 
boundary is shown in Figure (18) as an oblique solid line. As 
bar extension proceeds this boundary moves down the bar, closer 
to the center, material above the boundary recovering elastically, 
material below it continuing to deform plastically. The position 
of the boundary throughout the necking process qualitatively 
corroborates Bridgeman's [1] experimental observation that 
plastic deformation in necking bars is confined to the region 
between the inflection points of the neck boundary profile. 

The neck profiles shown in Figure (18) for necking in plane 
strain are in qualitative agreement with results obtained by 
Richmond £27] and Cowper and Onat [6] for rigid-perfectly plastic 
and hardening plastic materials, respectively. 

While the profile histories of Figure (18) provide a complete 
picture of neck development in plane strain this is not the case 
for plane stress. While the plane stress analysis proceeds in 
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the x-y plane of that figure it admits and indeed, in a thickness 
average sense, reflects the effect of out of plane deformation. 

The nature of this deformation is shown in Figure (19) where 
boundary profile histories in the three symmetry planes (x,y,z =0) 
are shown. It is apparent that the necking process in plane 
stress is properly considered as a problem in three spatial 
dimensions since a three dimensional neck develops. 

In plotting the profiles of Figure (19) the bar has been 
assumed to have an initial thickness of unity over its entire 
length. Thus the symmetry model for which results are shown 
is bounded by a symmetry plane at z = 0 and a traction free 
surface initially at z = 0.5. Note that the initial imperfection 
involved a reduction only in width at y = 0 not in thickness. 

The analysis does not predict the portion of the boundary 
profile at y = 0 (the root plane) which is nominally parallel 
to the z axis. The absence of shear in this plane and the 
symmetry condition at z = 0 suggest the profile shown. Merchant 
[31] has observed profiles of this nature in thin plate steel tensile 
specimens. 

The prediction of a three dimensional neck in plane stress is 
consistent with the more rapid load reduction and root plane stretch 
concentration previously noted. Figure (17), in this case as opposed 



Figure 19 Plane Stress Boundary Profiles 
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to that of plane strain. The result also exemplifies the limited 
utility of plane stress analysis which by definition neglects 
stress field three dimensionality which must follow from the 
deformed profiles of Figure (19). 

The Field Solution : The numerical results provide full 

histories of stress and deformation over the two dimensional 
domain of the analysis. The essential characteristics of these 
results are discussed below in terms of distributions of field 
quantities along the symmetry lines at x,y = 0. The numerical 
analysis provides stress component values at finite element 
centroids. The data plotted below represent averages of results 
for pairs of adjacent elements and are plotted at the centroid 
positions. No extrapolation has been performed. Thus results 
reported at x,y = 0 are actually values obtained at centroid 
positions slightly removed from these symmetry lines. 

Figures (20, 21) show histories of stress distributions at 
x = 0 in plane stress and plane strain, respectively. Distributions 
are given at maximum load, stage 1, as well as at stages 2, 3, and 
6 of the subsequent necking process. Associated x-y plane boundary 
profiles and elastic-plastic boundary locations are also provided 
for convenient reference. The plane stress and plane strain results 
are similar in form and variation during the necking process. 



Figure 20 Axial Stress Distribution: 
Plane Stress 
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Stress (lbs/in 2 xlO" 3 ) 




They are distinguished by the previously noted greater localization 
of the necked region in plane stress and the presence of a non-zero 
a z component in plane strain. A number of observations may be 
made which are equally applicable to both cases: 

1. The effect of the initial imperfection is apparent 
at maximum load, stage 1. Noticeable although small 
amplification of its effect upon the axial stress 
distribution has occurred. Recall that the initial 
local area reduction at y = 0 was 0.5 percent while 
at maximum load the loading direction normal stress 

0 ^ varies by 5 percent over the half length of the bar. 

The transverse stress a is non-zero only in the 

A 

vicinity of the initial flaw and attains a maximum 
value of approximately 2 percent of o^. 

2. At subsequent stages (2, 3, 6) of the necking process 
explicit correlation may be established between the 
axial stress distribution and neck geometry. All stress 
components respond dramatically to neck development. 

At stage 6 the root plane (y = 0) value of a (and in 
in plane strain a as well) is nearly twice that found 
at the opposite end of the bar. This distribution is 
directly related to the variation in cross-sectional 
area along the length of the bar. The root plane transverse 
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stress a , responding to boundary curvature in the necked 

A 

region, rises to 20 percent and 7 percent of a in plane 
stress and plane strain, respectively. Furthermore the 
sign of a in both cases correlates with the sign of 

A 

the boundary curvature, positive when the boundary is 
concave and negative when it is convex. 

3. The presence of elastic unloading behind (for greater 
y) the elastic-plastic boundary is clearly evidenced by 
the cusp in the o distribution occurring at the boundary 
location. The boundary position is also reflected in 
the relative magnitudes of the stress components Oy , a z 
at the various stages in the necking process. For 
example, in plane stress Oy at stage 6 exceeds at 
stage 3 in the plastic region ( y < 1.0) while the 
converse holds in the elastic region (y > 1.0). The 
progress of the elastic-plastic boundary through the 
bar may be followed by inspection of the axial stress 
distributions in the manner described above. 

The distribution of stress and deformation in the root plane 
of the neck is of particular interest since material property 
relations inferred from tensile data reflect average material 
behavior in this plane. Distributions of stress and loading 
direction stretch at several stages of necking are given in 
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Figures (22) for plane stress and (23) for plane strain. The 
field quantities are plotted as functions of x/x m where x m is 
the minimum width of the necked region. 

The distributions at maximum load, stage 1, again reflect 
the presence of the initial imperfection. In both cases variations 
in across the width is small (< 2 percent) while a x is no 
greater than 1 percent of at any point. As bar extension 
proceeds, however, the difference between necking in plane 
stress and plane strain is quite apparent. In plane stress the 
variation of stress and stretch across the root plane is 
significantly greater. 

The stress distribution in plane stress is the result of 
both diminishing thickness toward x/x m = 0, see Figure (19), 
and the presence of hydrostatic tension derived from the in 
plane boundary profile through the development of o x > 0 in 
the root plane. The a y variation is approximately 17 percent 
at stage 6. Note that the variation in stretch is noticeably 
less, about 10%, since it is the result of thickness variation 
only, plastic flow being independent of hydrostatic tension. 

In plane strain, on the other hand, variation of field 
quantities across the root plane results only from in plane 
neck geometry. This circumstance apparently provides sufficient 
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“stiffening" in the plane strain neck that boundary curvature 
remains much larger than in plane stress despite significant 
reduction in area at the root plane, 55 percent at stage 6. 
Variation in stress 5 percent, and stretch, 12 percent, at stage 
6 are much smaller than in plane stress. 

Figure (23) also provides a comparison of root plane stress 
distributions in plane strain predicted by equations due to 
Bridgeman [1] and by the present numerical analysis. Boundary 
profile data required as input to the Bridgeman analysis are 
extracted from the numerical results. Bridgeman' s equations 
are seen to overpredict a while underpredicting both a ^ and 
c x . While the absolute differences between the results are 
not large in the present case the trend of the comparison over 
the necking history considered suggests that for smaller root 
plane boundary profile radii the Bridgeman analysis may 
significantly underpredict hydrostatic tension.* The present 
results do not indicate under what circumstances such profiles 
might develop. 


*Needleman [30] draws a similar conclusion for the axisymmetric case. 



VI. INFERENCE OF STRESS-STRAIN RELATIONS 


The numerical data of Section V provide a basis for inference 
of stress-strain relations utilizing procedures similar to those 
employed for reduction of tensile data obtained experimentally. 
This exercise provides a basis for assessment of the validity 
of these relations since the actual stress-strain curve (111) 
for the material tested (analyzed) is known. Comparison of 
inferred and actual stress-strain behavior is undertaken over the 
full range of deformation from initial yield through development 
of a significant neck, stage 6 of the plane stress and plane 
strain necking processes of Section V. 

The comparison is performed in terms of effective stress 
and effective plastic strain. Procedures employed for the 
inference of these quantities from tensile data of Section V 
are described below. 

Vlam. S&i&te : In plane stress the effective stress is 

found as 



P /A 
y y 


where P is the axial load and A the corresponding minimum 

y y 

cross-sectional area of the tensile bar. The overscript bar 
in (116) denotes a cross-section average value. 


(116) 
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Prior to initiation of necking the average plastic strain in 
the bar may found in terms of overall bar extension as described 
in Appendix I. Thereafter, however, the overall extension is not 
indicative of minimum section, root plane, plastic strain. 
Alternatively, therefore, average minimum section plastic strain 
is computed directly from the minimum section area. The effective 
plastic strain F^ijPis given by 

7^)= ln\^= £n(A Q /A^ p )) 

in which incompressibility of plastic deformation has been noted. 

In (117) A is the undeformed area and A^ p ^is the deformed area 

° y 

corrected for elastic deformation. The corrected area A^ p )is that 

y 

which would exist should the applied load be removed and elastic 
recovery occur. In plane stress 

A^ p ) ~ A [1 + v 7/E] 2 

y y y 

where A is the deformed area under load. 

y 

It is apparent that for metals (a^/E « 1) the elastic area 
correction (118) will be small and of diminishing significance as 
plastic deformation increases. The data in hand confirm that for 
plastic strains in excess of several percent the elastic correction 
(118) is negligible thus permitting the use of A^t A^ in (117). 


(117) 
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Plane. S&uUn: Procedures similar to those described above 

are employed in plane strain. They are, however, complicated 
by the impossibility of inferring the plane strain (z) direction 
component of plastic strain from experimental data. Recall that 
in-plane strain 

d z = d^+ d(p)= 0 (119) 

where the elastic and plastic deformation rate components are 
not individually zero. To facilitate data reduction, however, 
we assume = /d z ^dt to be zero. The results below demonstrate 
this approximation to be of increasing accuracy as in- plane 
plastic deformation becomes large. The average effective plastic 
strain is thereby found as 

1.157 ln\f? )= 1.157 £n(A /A (p)) (120) 

ef y o y ' 

The elastically corrected area A ^ is found from 

A y (P) = Ay [ 1 + (v(l+v)/E) 0 y ] (121) 

which reflects the influence of the plane strain condition upon 
elastic recovery. As for plane stress the assumption A^t A 

y y 

is demonstrated to be a reasonable approximation for large 
plastic strains. 
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Computation of the effective stress is similarly impeded by 
lack of precise knowledge of a^. The ratio o z /o is given by 
Poisson's ratio v at the yield point and tends toward 0.5 for 
large plastic strains. Its variation with increasing deformation, 
however, cannot be ascertained experimentally. In lieu of such 
information we consider both limiting cases and find 


a e flv = 0,89 ° °y for CT z /a y = v = °* 3 

Oefl0.5 = °* 866 7 y f °r a z /a y = 0.5 

The present numerical results are within these bounds prior to 
necking initiation. 


( 122 ) 


Effective stress-plastic strain data inferred from the 
plane stress and plane strain results of Section V are given 
in Figure (24). The actual stress-strain relation (111) 
employed in the numerical analysis is also shown. 

Comparison of the inferred data points and the actual 
property curve for small plastic strains (<0.01 ) indicates, not 
unexpectedly, that noticeably more accurate prediction results when 
the data are corrected for elastic deformation. The corrected data 
for plane stress are quite accurate. Stress-strain data inferred 
from the plane strain results, however, are significantly in error 
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Figure 24 Inferred Stress-Strain Data 
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for effective plastic strains less than one percent. The error 
results from the assumption of zero plastic strain in the plane 
strain direction. This observation indicates that the inference 
of stress-strain data utilizing results of tensile testing of 
flat bars is not appropriate unless specimen dimensions ensure 
deformation under conditions closely approximating plane stress. 

For effective plastic strains exceeding two percent the 
results shown in Figure (24) clearly suggest that flat bar 
tensile data provide an adequate basis for inference of stress- 
strain data. The transition from the power law portion to the 
linear portion of the actual stress-strain curve, occurring 
for = 0.05, is accurately detected in both plane stress 
and plane strain. Furthermore it is evident that for plastic 
strains in excess of several percent the data need not be 
corrected for elastic deformation. The results also suggest 
that in the absence of necking the data reduction procedures 
employed are adequate for the inference of stress-strain data 
so long as the deformation is predominantly plastic. 

The effect of necking is evident in Figure (24). For 
effective plastic strains exceeding those corresponding to 
maximum load the effective stress is increasingly over-predicted 
in both plane stress and plane strain. At stage 6 of the necking 
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processes the error is approximately 10 percent in both cases, 
being slightly larger in plane stress. It is also clear, however, 
that the maximum load point does not provide a strict upper 
bound upon the utility of inferred stress-strain relations. This 
observation is consistent with the previously noted fact that 
necking does not initiate until somewhat after maximum load is 
attained. The data of Figure (24) indicate that reasonable 
prediction of material behavior may be obtained for effective 
plastic strains up to 50 percent higher than those existing at 
maximum load. 

The foregoing observations are, of course, strictly pertinent 
to testing of materials whose behavior may be characterized by 
stress-strain curves of the modified linear form (111). In 
particular it should be noted that in the above analyses the 
effective plastic modulus (da e j./d e ^^) is constant throughout 
the necking process. Consideration of a variable modulus 
would quantitatively affect the quality of stress-strain relations 
inferred from post-instability data. 

In order to demonstrate the effect of a ^ variable modulus 
the preceding plane stress and plane strain analyses have been 
repeated for the power law stress-strain curve of Figure (25). 

This curve is identical to the previously considered modified 
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linear relation for 0 < e('P)< 0.05. For e^> 0.05 the plastic 

ef ef 

modulus associated with the power law curve diminishes rapidly. 

Figure (25) provides a comparison of the two stress-strain curves 
as well as their associated plastic moduli. 

Effective stress-strain data inferred from the post- instability 
results of plane stress and plane strain analyses are compared 
with the power law input relation in Figure (26). Two significant 
effects of the variable and diminishing modulus are immediately 
apparent. Maximum load is attained at an effective plastic strain 
much lower than that found in the case of a constant plastic 
modulus. Furthermore* the error in stress-strain curve prediction 
subsequent to attainment of maximum load develops more rapidly 
for the power law material than is observed in Figure (24) for 
a linear hardening material. 

The observed effects of a variable plastic modulus are 
qualitatively consistent with the nature of the equations governing 
the finite deformation process. Inspection of the velocity 
equilibrium equations (67) reveals that the nature of the deformation 
process is controlled by the relative magnitude of the existing 
stresses and the material stiffness. In the present case of a 
continuously diminishing modulus the governing equations are 
dominated by the effect of existing stress for smaller total 
deformation than would be the case for a constant modulus. Thus 
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maximum load is attained for smaller uniform bar extension as 
observed in the numerical results and predicted by the Considere 
criterion (113). Subsequent to necking initiation the plastic 
modulus of material in the necking region diminishes with increasing 
overall bar extension. Thus a more localized neck, and consequently 
higher root plane hydrostatic tension, will develop accounting for 
a larger error in inferred effective stress for the power law 
material . 

The preceding results suggest that stress-strain relations 
inferred from flat bar tensile data are highly accurate only for 
a bounded, material dependent range of plastic strain. This 
range excludes both the vicinity of initial yield, wherein 
non-proportional loading renders the inferred effective quantities 
indeterminate, as well as the large strains associated with tensile 
necking. 

Of equal significance, however, is the demonstration of the 
utility of finite deformation analysis capability in assessing the 
validity of stress-strain relations developed from test data. The 
present solution capability admits consideration not only of 
functional relations such as those considered above but also of 
numerical relations provided directly by experimental data. 
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The analysis might profitably be employed for evaluation of more 
sophisticated property inference procedures as well as for the 
prediction of experimental load-deformation data on the basis of 
inferred stress-strain relations. The latter approach provides 
a direct means of assessing the validity of stress-strain relations 
subsequently to be employed in design analysis. 
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VII. CONCLUDING REMARKS 


A theoretical basis has been established for analysis of 
finite deformation of metals. The observation that finite 
deformation of such elasto-plastic materials may be viewed as 
a process rather than an event has led to the derivation of a 
complete initial- and boundary-value problem distinguished by 
its quasi -linear nature. This feature of the formulation motivates 
the adoption of an incremental approach to numerical problem 
solving. 

Efficient numerical solution capability has been developed 
for problems of two dimensional deformation under conditions of 
either plane stress or plane strain. The validity of the 
numerical analysis has been evaluated by considering a variety 
of elastic and elasto-plastic finite deformation problems whose 
homogeneous nature renders analytic solution possible. It is 
demonstrated that accurate solutions may be obtained for problems 
involving extremely large displacements and rotations. 

The numerical analysis has been employed for the investigation 
of necking in flat metal tensile bars. The results of this 
investigation provide not only the first full numerical solutions 
for tensile necking of metals in plane stress and plane strain 
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but also an appraisal of the validity of stress-strain relations 
inferred from tensile test data. It is demonstrated that such 
relations inferred from the behavior of flat tensile bars are 
erroneous for both very small and very large (post-instability) 
plastic strains. The post-instability error is shown to be 
significantly dependent upon material behavior and in particular 
upon variation of the plastic modulus. 

It is evident from the results obtained that present knowledge 
of the mechanics of tensile testing is insufficient to enable 
precise characterization of material behavior from tensile data 
over the full range of a test. The theory and numerical analysis 
which have been developed- provide the means for necessary further 
study of tensile testing mechanics and procedures. Such 
investigations might consider, for example, the effects of material 
properties and tensile bar geometry upon the necking process. 

The possibility exists of developing procedures for correcting tensile 
test data to account for root plane hydrostatic tension and thereby 
to provide a basis for inference of accurate stress-strain relations 
from post- instability test results. The existence of such a 
material independent correction procedure is suggested by the 
experimental work of Bridgeman [1]. 

The availability of analysis capability for finite elasto-plastic 
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deformation also provides the basis for evaluation of alternative 
material property tests. Two obvious candidates present themselves, 
compression testing and indentation testing. Both of these are 
attractive since neither involves an instability phenomenon. Since 
full solutions for stress and deformation may be obtained the 
inhomogeneous character of the indentation or hardness test would 
not necessarily obstruct inference of effective stress-strain 
data. The utility of these alternative procedures might be assessed 
not only by analysis of the individual tests but also by comparison 
of stress-strain data inferred from the results of analysis of 
several test methods. 

The present analysis of necking in flat tensile bars, as well 
as the possible avenues of research identified above, suggest the 
primary significance of the finite deformation solution capability 
which has been developed. Precise characterization of inelastic 
material behavior can be extracted from the results of mechanical 
testing only if the mechanics of each test employed is understood 
and, in particular, only if quantitative distinction can be made 
between the effects of material and geometric nonlinearity. The 
availability of an analysis incorporating treatment of both forms 
of nonlinearity provides the means of distinguishing their effects. 
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APPENDIX I 


Operational Form of the Elasto-Plastic Constitutive Equations 

In order to utilize the elasto-plastic flow equations (50) 
and the inverse equations (53) we must choose a specific form 

for the yield function <f> = t of (39) and operationally define 

the equivalent plastic modulus y ^ of (46). Means must be 

provided for the evaluation of y (p) from test data for 

particular materials. 

In the present analysis <#> is taken as the octahedral shear 
stress r 

o 

♦ - T - T 3 [(2/3) J,] 1/2 
eq o 2 

where J. is the second invariant of the deviatoric stresses s. 

2 ij 

J 2 - 0/2) sU Sij 

An equivalent plastic strain rate d eq (P) is defined (43) by 
requiring the rate of plastic work W^'to be given as 

w(P) = t d 

o GC| 

from which we find 

d (P) = 3d (p) 

eq o 


(I-D 


( 1 - 2 ) 


(1-3) 


(1-4) 
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In (1-4) d (p) is the rate of octahedral plastic shear strain 

o 

d (p) = [(1/3) d..(P) d i J(P)] 1/2 (1-5) 

o ij 

Integration of (1-4) with respect to time defines the equivalent 

plastic strain e (p) in terms of the octahedral plastic strain 
eq 

Y (P). 

O 

Y Efd^dt 

° J ° (1-6) 

e ^ = f d ^ dt = 3 /* d (p) dt = 3 y (p) 
eq J eq J 0 'o 

The equivalent plastic modulus y g q^ P ^ ma y now be expressed in 
terms of an octahedral plastic shear modulus y 

o 

y eq (P) = (1/2) dx eq /d Eeq (P) = (1/3) y# (P) (1-7) 

y (p) = (1/2) dt /d Y ^P^ (1-8) 

o o o 

The octahedral plastic shear modulus may be evaluated for 
particular materials utilizing data obtained by quasi -static 
testing under simple loading conditions. In the case of uniaxial 
tension* we have a single non-zero stress component a Y (t) > 0 
and from (1-1) we find the octahedral shear stress 

x = (.£73) o x (1-9) 

O A 


*The discussion is restricted to tensile data obtained prior to 
initiation of necking. 
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Note that a is defined as applied load divided by actual specimen 

A 

cross-sectional area as discussed in detail by MacGregor [32]. 

As a x increases with time the tensile specimen deforms in simple 
extension. For an isotropic material this homogeneous elasto- 
plastic flow process is described by the cartesian deformation 
rate components 

d = d (p) + d ( e ) > 0 
xxx 

d = d (P) + d ( e ) < 0 

y y y (1-10) 

d = d (P) + dJ e ) < 0 

Z 2 z 

d = d 

y z 

Noting the incompressibility of plastic flow we find 

dy (p) = d z (p) = -(1/2) d x (P) (I-ID 

The plastic strain rate in the loading direction is given by 

d^(P) = U/£) - a/E (1-12) 

where l is the instantaneous x dimension of the tensile specimen 
gage length and E is the elastic Young's modulus. Substituting 
(1-12) into (1-5) and integrating with respect to time the octahedral 
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plastic strain is found as 

Y (P) - (1/^) IU\ - a /E] (1-13) 

o X * 

where A = l/l and l is the original (undeformed) gage length.* 

X o o 

The octahedral property relation x (y ^ ) provided by 

o o 

(1-9) and (1-13) may be approximated by a monotonic function or 
retained in the numerical form provided by the experimental data. 

The FIPDEF program utilizes the octahedral data directly and 
employs a finite difference technique to re-evaluate y in 

O 

each finite element at the beginning of each time step. The 
value of x for which y ^ tends to zero defines the initial 

O o 

proportional limit for the material. Elastic analysis is 
performed until the octahedral stress exceeds this value. 

It is occasionally convenient to utilize material property 
variables which reduce to the principal uniaxial quantities in 
the case of simple extension under uniaxial loading. For this 
purpose we define 

°ef = [3 J 2 ]!/ 2 

d ef ( p > = [(2/:3) d f >) (1-14) 

V (P> ■ / d ef (p)J<it 

♦Under the assumption of infinitesimal elastic strains («1) the 
undeformed length defines an appropriate reference state for the 
entire plastic deformation process. 
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In the case of uniaxial tension these reduce to 


a x = a 
ef x 


d (p) = d (p) 
ef x 


(1-15) 


(1-16) 


These effective quantities (1-13) are related to the octahedral 
variables as 

t = (/2/3) a x 
o ef 

d (p) = (1//2) d (P) 
o et 

= (1/3) da e f/ de e f^ 

The elasto-plastic flow equations (50) and inverse equations 
(53) may be written in terms of the octahedral quantities as 


2jidi = o\ -(y/l+v) o\ 6i +(v/p (P>)(3x 2 )' 1 s i . s kl c, , (1-17) 

J J K J o o J Kl 

a 1 '-? = Xd k . 6 1 ' .+2pd 1 ' .-2p[3x 2 (l+u ^ p Vp)] ^ s 1 . s k ^ d (1-18) 

JPk Jjoo jkl 


where x,y are the Lame constants of linear elasticity. 
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APPENDIX II 


Elasto-Plastic Constitutive Equations for* 
Plane Stress and Plane Strain 


The general elasto-plastic constitutive equations of Section 
1.2 are specialized in Appendix I for analysis of materials in 
which plastic yielding is controlled by the octahedral stress. 
These equations (1-17, 18) are expanded below for analysis of 
elasto-plastic flow under conditions of either plane stress 
or plane strain. 

It is convenient to develop the equations in matrix form. 

For this purpose we define matrix vectors consisting of the 
in plane cartesian components of the Jaumann stress rate a.. 

‘ J 

and deformation rate d^j. 


a 


C 



0 



(II-l) 


(II-2) 
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For either plane stress or plane strain the elasto-plastic flow 
equations may be written 


DS = B W.d 5 ; 5, C = 1, 2, 3 ( I 1-3) 

and the inverse equations as 

a 5 np??D5 ; t, 5*1,2,3 (II-4) 

The constitutive matrices and P ? 5 are full and symmetric. 

Specific forms of these matrices are developed below. The 
notation employed is similar to that utilized by Swedlow C 33 ]- 

Plane Stress : It is assumed that 


a = a = a =0 

z xz yz 


( I I -5 ) 


Under these conditions the constitutive matrices are found as 


= O/E) 


T + s x 2 /s 2 

X o 


(sym) 


- v + s xV s . 2 
1 + s y 2/s . 2 


2s a /s 2 
x xy o 

2s a / s 2 
y xy o 

2(l+v+2a 2 /s 2 ) 

xy o 


( I I— 6 ) 
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pec = E/(l -v 2 ) 

1-P(s +vs ) 2 /s 2 v-P(s +vs )(s +vs )/s 2 (v-l)P(s +vs )a /s 2 
xy 0 xyyXo x y *y ° 

l-p(s +VS ) 2 /S 2 (v-l)P(s +vs )cr /S 2 

y x o y x xy o 

(sym) [(1-v)/2]-(l-v 2 )Pa /s 2 

xy o 

where 

s 2 = 6 T 2 v ( P)/E ( I 1-8 ) 

o o o 

1/P = l-v 2 +[s x 2 +2vs x s y +s y 2 + 2{l-v) Cxy 2 ]/S o 2 ( I I -9 ) 

In plane stress the octahedral stress is found from (1-1) and 
u ^ is the octahedral plastic shear modulus (1-8). 

O 

t o 2 = (2/9) (o x 2 - a x o y + o y 2 + 3 (11-10) 

The deviatoric stresses are 

s * = (2o x ■ V ' 3 

s = (2c„ - c ) / 3 (11-11) 

y y x 

S 2 * ‘ (c X + °y> ' 3 


(11-7) 
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The constitutive equations (1-17) also provide the auxiliary 
relation 

d z - (1 /e) t(-v+s z s x /S o 2 )5 x + (-v+s 2 s y /s o 2 )y 2(s z c xy /s o 2 ; xy )] 

Plane Strain : It is assumed that z/dz is a null operator and 

that 



Under these conditions the constitutive matrices become 


E(l+s 2/s 2) 

Z ° 

1-v 2 +(s x 2 +2vs s +s 2 )s 2 -v(l+v)+(s s -2vs 2 )/s 2 2[(s +vs )a /s 

A x^zo x y Z O x z xy 

l-v 2 +(Sy 2 +2vs s +s 2 )/s 2 2 [ (s +vs )a /s 

y y Z L o yZ^y 

(sym) 2{(l+v)+2[a xy 2 +(l+v)s z 2 /s 

P ?? = l/(l+v) 


3k(1-v) - Q s 2 /s 2 

X O 

3kv - 

Q VA* 

- Q S xV s 


3k(1-v) - 

Q s y 2 /S Q 2 

^ s y a xy /a 




-Q a 2 /s 2 




Xy o 


( 11 - 12 ) 


(H-13) 


2 

2 (11-14) 


(11-15) 


144 



where s 2 is defined by ( I 1-8) and 


K " °y " a z } ' 3 


In (11-15) 


s „ “ (2a 


°z " a x } ' 3 


2 = 


(2CT Z " °X " CT y } / 3 

(2/9)(a x 2 +a y 2 +a z 2 -a y a z -a x a z -a x a y + 3a xy 2 ) 


(11-16) 


Q = 2p/(1 + u /p (p)) 


(H-17) 


and k is the elastic bulk modulus. 


k = E/ 3(1 -2v) (11-18) 

Note that boundedness of k and thereby P ^ in (11-15) requires 
v < 1/2. Hence plane strain analysis of elastically incompres- 
sible materials is not possible. Approximate analysis 
of such materials may be accomplished by setting k/E »1 . 

In addition to (11-14, 15) we have the auxiliary relation 


o = a - 1/(1+v){[3kv -Qs S /s 2 ]d + [3kv-Qs s /s ^]d 

z z 2 X o X 2 y o y 

- [2Q S a /s 2 l d > 
z xy o X y 


(H-19) 
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APPENDIX III 


The Rate Stiffness Matrix for Planar Analysis 

The rate stiffness equations are specialized below for analysis 
of planar problems under conditions of either plane stress or 
plane strain. A triangular finite element is employed with the 
assumption that the velocity field varies linearly within the 
element. 

The Linear Velocity Element : Consider the triangular finite element 

of Figure (III-l). The instantaneous configuration of the element 
is defined by the coordinates of its vertices or nodes. These 
nodal coordinates are represented in matrix form as X a 



where X^ (Y^) is the x (y) coordinate of node 1, etc. 



Figure III-l Linear Velocity Finite Element 
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Vectors of nodal velocity components V a and nodal loading 
rate components T a are similarly defined. 




148 



The velocity field components within the element are 
approximated as 

v x = A“ I x a (x,y) 

a 1 5 • • • 56 

Vy = A a * y a (x,y) 

where the vector functions <f> a possess the components 


a 

1 

2 

3 

4 

5 

6 


1 

0 

x 

0 

y 

0 


t a 

y 

0 

1 

0 

x 

0 

y 


Hence the velocity field approximation ( I I 1-4) is expanded as 

\ (x,y) = A-| + A 3 x + A 5 y 
v y (x,y) = A 2 + A 4 x + A g y 

Evaluation of the velocity field approximation ( 1 1 1 -6 ) at the 
nodal positions (III-l) must yield the nodal velocities ( I I 1-2) 


(II 1—4) 


(III -5 ) 


( 1 1 1 -6 ) 
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Substituting (III-l) into (III-6) and solving the resulting 
equations for A 01 obtains 


A a = r a6 V e a, 3 = 1 ,...,6 ( I I 1-7 ) 


r a& = (1/2A) 

( X 2 Y 3' Y 2 X 3 } 0 

o (x 2 y 3 -y 2 x 3 ) 

0 


(V y 3 ) 

0 


( Y 2 - Y 3 ) 


(X 3 -X ) 0 


( x 3 y 1 " y 3 x 1 > 0 (XiY 2 -Y 1 X 2 ) 0 

0 ( X 3 Y r Y 3 X l) o (x 1 Y 2 -Y 1 x 2 ) 


( Y 3 " Y 1 ) 
0 


(x r x 3 } 


( x 3 - x 2 ) 


and A is the area of the element. 


0 

(VV 

0 

( x r x 3 ) 


( y i- y 2 ) 

0 

(X 2 -X 1 ) 

0 


0 

< Y T Y 2> 

0 


(x 2 -x 1 ) 


A = (1/2) [X 2 (Y 3 -Y 1 ) + X-| ( Y 2 “ Y 3 ) + X 3 (Y 1 -Y 2 )] (III-9) 


Substituting ( 1 1 1-7 ) into (III -4) yields the velocity field 
representation ( I II - 1 0 ) corresponding to the general form ( 74 ) 
of Section III .1 . 

The element deformation rate and Jaumann stress rate fields 
in the x-y plane are given by the matrix vectors and D e , 
respectively (? = 1,...,3) 


(III - 8 ) 
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(III-10) 


A 

a 


d 

X 


X 

A 

. a 

D ? = 

d 

y 


y 

A 



a 


2d vw 

_ xy 


xy 


These in-plane rates are related by the constitutive equations 
of Appendix II 

c ? = D 5 ; a, g = 1,...,3 (III-ll) 

where P ^ = P^ ? takes a form appropriate for either plane 
stress ( I 1-7) or plane strain (11-15). 

The Planar Rate Stiffness Matrix : The preceding formulation of 
the rate behavior of the linear velocity element provides the 
basis for specialization of the element rate stiffness matrix 
K a & of (79) to a form appropriate for two dimensional analysis. 

The general form of this matrix is given below employing cartesian 
tensor notation*. The matrix is decomposed into three components 
to facilitate the ensuing two dimensional specialization. 

*The notation of Section III is employed here. The distinction 
between covariant and contravariant tensor components need 
not be made in cartesian coordinates. 
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(III-12) 


K a6 




K “S 
2 



+ i^ae + K aB a, g = 1 6 

/ "“S, p ijki *ki ^ 



T 

{r aT1 [(j) n 




a. 9. 

IP J»P 



B 

m 



a. <|>. 6 .r 6e }dl/ 

IP J.J 


jm 



]r 6 e}dl/ 


In (III-12) matrices and are symmetric while K.j afJ is not. 
Latin indices take the values 1, 2, 3. 

Expansion of (III-12) for analysis of plane stress and plane 
strain deformation procedes on the basis of the following 
observations. 

1. Since all quantities present in (I 11-12) are uniform in a 
single element the requisite volume integration over B m reduces 
to multiplication by 1/ the volume of the element. 


‘'n, ' V 


( ni - 13 ) 


In ( I I I -1 3) A is given by (III-9) and l , the element thickness 

m 

in the z direction is computed by noting that in plane stress 


d 

z 


= 1/1 
m m 


(III-14) 
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(III-15) 


where is given by (11-12), while in plane strain 


An 


= 1 


2. Substituting for in the matrix may be written as 


V s f T 

' J P^ f ^ 6 


dl / 


oc , |3 1 , . . . ,6 

5 » ? = 1> 2 , 3 


m 


where P ^ is the constitutive matrix in (III-ll) and 


T° 5 =(1/2A)- 


( Y 2-Y 3 ) 0 (Yg-Y-,) 0 (Y r Y 2 ) 0 

0 (X 3 -X 2 ) 0 (X-, -X 3 ) 0 (X 2 -X 1 ) 

[(X 3 -X 2 ) (Y 2 -Y 3 ) (X r X 3 ) (Y 3 -Y 1 ) (X 2 -X-,) (Y r Y 2 )J 


3. Reduction of to proper form for planar analysis is 
accomplished by considering Latin indices to take the values 
1,2, corresponding to the x and y coordinate directions. 

4. The matrix contains a factor <f> . . associated with 


J,J 


the dilatation rate v. ., i.e., 

J 5 J 


(III-16) 


(III-17) 


. r 6a V a =v--=d +d + d 

9 j,j J.J x y 2 


(III-18) 
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In plane strain v. • is completely defined by the in- plane 

J 9 J 

deformation rate components d and d (d =0) and K, a 0 may be 

x y z J 

employed as written in (III-12) for i,j,p = 1,2. In plane 
stress, however, d z / 0 and is given by (11-12). To accommodate 
both of the above cases the matrix is written as 


where 


1/ a& 

S 


■/ 


I r a n ^ n p* r «e}dt/ 

i >P 


m 


P 6 = (l/3ic) 


E 11 + E 13 


E 13 + E 23 


E 13 + E 23 


E 12 + E 22 


For analysis of plane stress E in ( 1 11-20) are the P^ ? of 
the constitutive equations (III-ll); for plane strain we take 


E 11 = E 22 = 3k 

E 12 T E 13 ‘ E 23_ = J?-. - 


C 1 11-19)- 


( I 11-20 ) 
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and in both cases k is the elastic bulk modulus k = E/3(l-2v). 


In summary the rate stiffness matrix for planar problems is 
written as 


where 


K“ e = + k 2 “ 3 + K 3 a3 


(II 1-21 ) 


ag 


= (£ m A) {P 5 Pte r ?6 } 


ag 


K 3 a3 


= U A) {rn* 1 ? . a, „ -26?. a. f? ]r^> 


m 


j,i ipJ.P 
= (1 A) {r“V P 6 r 6e } 

m i >p 


ij jm im 


In ( I 11-20) Latin indices take the values 1, 2 while of the 

Greek indices; take the values 1, 2, 3, while the remainder 

vary over 1,...,6. Analysis of either plane stress or plane 

strain procedes by choosing appropriate forms for P ^ and P 5 

and setting the value of l according to either of ( II 1-1 4) or 

m 

(III-15). 
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APPENDIX IV 


Solution of Verification Problems 

Homogmejcnu, FZnvte. V&^oumcution 

Problems of homogeneous finite deformation may be defined 
by specification of velocity fields v. (x-l.t) whose gradients 
are spatially uniform, that is, 

(IV 

The time dependent stress fields corresponding to velocity 
fields of the form (IV-1) are likewise spatially uniform. 
Consequently the velocity equilibrium equations (67) are 
satisfied identically by any velocity field corresponding to 
homogeneous deformation. 

Complete solutions to such problems are developed by time 
integration of the velocity field to define the deformed 
configuration of a body and of the constitutive equations (53) 
to define the stress field in the body. Solutions are developed 
below for problems of finite homogeneous extension and simple 
shear. Except as noted the solutions are developed in a fixed 
cartesian coordinate system (x,y,z). 



A-IV.l Finite Extension 


Consider the problem of extension of an isotropic material 
in the x-y plane as defined by the cartesian velocity components 


v x = ox 
V y = 


(IV-2) 


where a is a constant and k is an unknown function of the material 
properties and possibly the state of stress. We may consider 
cases of plane strain for which 


v z - 0 


(IV-3) 


and of plane stress for which we expect 


v z = kaz 


( IV-4) 


The deformation rate components for this velocity field are 


d “ a 
X 




d z ’ 


0 plane strain 

ka plane stress 


d xy d xz d yz = 0 


( IV— 5) 
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Integrating (IV-5) with respect to time the total deformation may 
be described by the logarithmic strain components 

e = £n. X x = at , 

e. f = In x w = Un X (IV-6) 

y y x 

0 plane strain 

e = In X = 

1 K£n A plane stress 

A 

In (IV-6) (X Y ,x ,x ) are coordinate direction stretch ratios 

and K is expected to be a function of material properties and 

possibly the state of stress. 

Non-zero stress components a , a and in plane strain a 

A y z 

are found by substituting the deformation rate (IV-5) into 
the constitutive equations (53) and integrating with respect 
to time. The Jaumann stress rate. a., in (53) reduces to a 

* J 

material derivative since no rotation occurs. 

In plane applied loads are found by integrating the tractions 
t Y and t over x and y coordinate faces of the body, respectively. 

Ay 

That is, for a body which is initially a unit cube, as in 
Section IV, we find 


P = a A on x faces 

X XX 

P = a A on y faces 

y y y 


(IV-7) 
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where 


A = X X 
x y z 


A 

y 


X x X z 


Elastic Material : Substituting the deformation rate components 

(IV-5) into the constitutive equations (53) for v>/u e q^ = 0 
yields 


° x = U+2y) d x + X(dy+d z ) 

o y = U+2y) dy + X(d x +d z ) 

° z = U + 2y) d z + X ( d x +d y) 



Integrating (IV-9) with respect to time for an initially 
stress free material obtains 


a x = ( x+2u ) e x + x ^ e y +e z^ 

a = (X+2y) e + X(e +e z ) 

y y x Z 

o = (X+2y) e + X(e +e ) 
z z x y 


where the strains {e , 

x 



e z ) are given by (IV-6). 


( IV-8) 


(IV-9) 


(IV-10) 
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Solutions to problems of plane stress and plane strain extension 
may be given in terms of a prescribed stretch ratio in the x 
direction X = X. 

A 

Plant S&LZA& 

a - E In X 

A 


Plant SlAaln 



0=0 

z 


X = x~ v 
z 


X X = X" 2v 

y 2 


XX- X 1 “ v 
X Z 


o x A x * E X" 2v InX 


\V° 


<J = [E/ {T-v 2 ) ] £n X 

A 

a z = [ E/(l~v 2 )] In X 
°y = 0 

X y = xCv/ M)] 

X = 1 


(iv-11) 


(IV-12) 
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A 

y 


p 

X 


p 

y 


yz = x ~[v/(v-D] 

A x A z = A 

a x A x = [E/(l-v 2 )] ^/(v-l)]** X 


Elasto-Plastic Material : Development of solutions for elasto- 

plastic materials is facilitated by considering the flow 
equations (51) rather than the inverse equations (53). 
Substituting the deformation rates (IV-5) into these equations 
and noting that the shear stresses are zero yields 


d x = ( 1//E ^ a x" v ^ a y +a z^ + ^ 3>/2 ^ a ef' /3a ef^ a x"^ y/3 ^ a x +a y +0 z^ 

d = (1/E)[a -v( a +a ) ]+ ( 3/2 ) ( a f / 3a f ) [ct — (1/3) ( a +o +a ) 
y y y ^ eiciy x y ^ 

d z = (l/E)[a 2 -v(a y +a z )]+(3/2)(a ef /3a ef )[a 2 -(l/3)(a x +CT y +a z ) 
d xz = V = d xy = 0 

where is defined by i-l 4 and in the present case is given as 
o ef 2 = (l/2)[(a x -Oy) 2 + (° y -° z ) 2 + (a z -a x ) 2 ] 


(IV-12) 


(IV-13) 
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In (IV-13) we have 


1/3= 0 for < ay ; elastic flow 

3 = da e f/ cle e f^ P ^ > 0 for a e f - a V > elasto-plastic flow 

where Oy is the initial yield stress of the material. 

Immediate integration of (IV-13) is not possible due to 
the presence of the stress components in the equations. To 
facilitate integration of these equations we assume proportional 
loading, i.e.. 


a ij = q(t) a ij° (IV-14) 

where a .- 0 is a constant reference stress state defining stress 

I J 

component proportionality. The time varying loading level is 

controlled by the scalar q ( t) . Integration of (IV-13) subject 

to CIV-14) provides equations relating total stresses and 

logarithmic strains (IV-6). This simplified formulation, termed 

deformation theory, is correct to the extent that the constraint 

(IV-.14) is valid for particular problems. For plane stress 

extension the stresses have the form (IV-14). In plane strain, 

however, the proportionality between the in plane stress 

components (a ,a ) and out of plane component a changes at the 
x y z 
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yield point. Hence the deformation theory is in a strict sense 
not appropriate for plane strain elasto-plastic analysis. For 
the extension problem considered here the quantitative error 
in deformation theory results is small and the availability of 
those results in analytic form facilitates comparison with 
FIPDEF numerical results. For more general problems, however, 
the general flow theory of Section 1.2 must be employed. 

For < ay the solutions for plane stress and plane 
strain are given by (IV-11) and (IV-12), respectively. In 
each case a yield point value of X, x=Xy, may be found 
corresponding to a =Oy. For > ay the solutions are as 
follows. 


Plane. S&uu>t> 

X > Xy = eYTi 

CT V ■ [3/(l+n)][£nX+y] 

J\ 

°y = <, z = 0 

P = {[B/(l+n)UttX~+Y]/X> e r o 
x 

P = 0 

y 

where e is the natural base. 


( IV-15) 


163 



In ( IV- 15) 


n = 3/e 
Y = CTy-/ 3 

To = [n(l-2v) (£nX+y)]/(l+n) 

Plane. Strain 

X > X y = e ^ n(1_v2 ^ 

In X = H-. a v + H 9 
y • x 2 

a x = Hg In X + Hy 

”y =o 

a z = CT x n+2vn)/2(l+n)] + cyf>[(2v-l )/2(l+n)] 

P x ■ H 3 x y ** x + H 4 x y 

where e, n» y are as defined previously and 

H.j = -(l+2vn) [3+2n(l+v)] / [4(l+n)3 ] 

Hg = [y<I , /4(1+ti)] [4vn+2n+3+4nv 2 ] 

H 3 = 43(l+n) / [3+4n(2-v)+4n 2 (l-v 2 )] 

H 4 = (4n+3<(>^) / [3+4n(2-v) + 4n 2 (1-v 2 )] 

l/4> 2 = 1-v + v 2 


(IV-16) 
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Maxi mum Load Phenomenon : The. applied load in the direction of 


extension may be written as 

P = a A = f (X) (IV-17) 

x x x 

This relationship P (X) is characterized in all cases considered 

A 

here by a maximum value P £ occurring at some critical stretch 
x c . The critical stretch may be found by solving 

d/dx (P Y ) = 0 (IV-18) 

A 

where P (a) is given by any of (IV-ll, 12, 15, 16). The maximum 

X 

load is evaluated by substituting x = X in the appropriate 

0 

load equation. 

Expressions for P , x are given in Table (IV-1) for each 
c c 

of the four extension problems solved in this section. 
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TABLE (IV-1) 
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Young's Modulus n = 3 Ti = l+nv 2 -(l+2vn)^/4(l+n) 

Poisson's Ratio Y = <We 

da-/de e (p) = constant r = (1+n)/(l+2vn) r 2 = [3+4n(2-v)+4n^(l-v )]/{(l+2vn)[3+2n(l+v)]} 

uniaxial yield stress 1 /<j> 2 = 1-v+v 2 

natural base r 3 = (4n+3<t> 2 )/[4(l+n)] 



A-IV.2 Unidirectional Extension 


Consider the homogeneous deformation corresponding to the 
velocity field components 


v x i = ax' 



( I V - 1 9 ) 


where (x 1 , y', z') is the rotating coordinate system of Figure 
(10b). The deformation rate components in the rotating system are 


d x , = a 

dy* ^z 1 ^x'y 1 ^x'z' ^y'z 1 ® 

and the total deformation is described by the single non-zero 
strain component 


(IV-20) 


e x . » In \ x , (IV-21 ) 

For an elastic material the stress field in the prime 
system is found by writing (IV-10) in that system and setting 
all strains to zero except e , . The complete solution of the 

A 

problem in (x', y', z * ) is 
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X = X 

x 1 



cr v i = { (1-v) E/[(l+v) (l-2v)]} In X 
x (IV-22) 

o y i = a z i = (v/l-v) a x 


P x . = V A x' = VVV = 1 


Pyi ~ a y'Ay' = OyiX “ (V/1-V) X a x i 


where the initial configuration has been taken as a stress free 
unit cube. 
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A- IV. 3 Simple Shear 


The velocity field (IV-23) is imposed upon an initially 
stress free elastic continuum. 


v x = 2ky 



k = constant 


(IV-23) 


Corresponding to (IV-23) are the non-zero deformation rate 
components 



and non-zero spin components 


( IV-24) 


“xy = -£0 yx = k ( IV-25) 

Recalling the definition of the Jaumann stress rate, a ^ 

5 k = CT k + a k a ) m . - a m , io k (IV-26) 

1 l m 1 i m 

The constitutive equations, (53) for y/y ^ = 0, may be 

eq 

wri tten 
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(IV-27) 


d/dx(cr ) = a 

X Xy 

d/dx(a y ) = >a xy 
d/dr(a xy ) = U + (1/2) (a y - a x ) 

where : 

2y = E/(l+v) 

x = 2kt 

Integrating (IV-27) subject to initial conditions of zero stress 
obtains the solution 

o = u sin t 
xy 

a x = y (1-COS x) (IV-28) 

a y = y (COS x - 1 ) 
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